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I. INTRODUCTION 



A. DESCRIPTION OF THE PROBLEM 

A statistical design is a plan according to which an experiment is patterned. It 
provides the basis upon which appropriate statistical tests and inferences can be made 
after the experiment has been performed. The selection of the experimental design to 
be used in a given situation is extremely important because it plays a predominant role 
in the efficiency of the experiment, the precision with which the objectives are met, and 
the total effort (and cost) expended upon the experiment. 

The concern of this thesis is the development of power algorithms for several 
common tests, including the F-test and t-test. The experimenter wants a test to 
achieve the correct decision with as high probability as possible in practical operational 
tests and evaluations given certain conditions hold. 

B. SCOPE OF THE THESIS 

An interactive computer program is developed which uses power algorithms for 
several common tests. The program can be used by decision makers who may not 
have deep knowledge of statistics or may be unfamiliar with the details of the design of 
experiments. The output helps the decision maker design experiments appropriate for 
testing hypotheses with common statistical tests. The experiment designer can use the 
program to evaluate the effect upon power of varying underlying design parameters 
such as sample size. 

Because the underlying statistical noncentral distributions used in computing 
power do not have tables, approximation methods are used for computing the 
noncentral CDFs. Approximation methods were developed which are efficient, 
requiring reduced computer CPU time. 

C. BACKGROUND 

A hypothesis is a statement about the values of the parameters of a probability 
distribution. For example, suppose we think that the mean yield of a chemical process 
is more than 94.5 percent. Hypotheses to test this statement might be expressed 
formally as 

Ho : H < 94.5 
Ha ; n > 94.5 
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The statement Ho : ^ 94.5 is called the null hypothesis, and Ha : ^ > 94.5 is called 

the alternative hypothesis. The value of the mean specified in the null hypothesis 
might typically be determined in one of two ways. It may be the result of some theorv' 
or model regarding the process under study, or it may be the result of contractual 
specifications. 

To test a hypothesis one usually devises a procedure for taking a random sample, 
computes an appropriate test statistic, and then rejects or declines to reject the null 
hypothesis Ho, depending on the outcome on the test statistic. Part of this procedure 
is the specification of the set of values for the test statistic which lead to rejection of 
Ho. This set of values is called the critical region or rejection region for the test. 

Two kinds of errors may be committed when testing hypotheses. If the null 
hypothesis is rejected when it is true, then a type I error has occurred. If the null 
hypothesis is not rejected when it is false, then a type II error has been made. The 
probabilities of these two errors are given special symbols: 
a = P(type 1 error) = P(reject Ho | Ho is true) 
p = P(type II error) = P(fail to reject Ho | Ho is false) 

Generally, the determination of p requires the use of a noncentral sampling 
distribution (e.g., noncentral F, noncentral x ^ , ^ind noncentral t) for which tables are 
not readily available. The noncentral sampling distribution depends on a parameter 
called the noncentrality parameter. The noncentrality parameter is usually a measure 
of the distance (in some sense) between the values of the parameter under the null and 
alternate hypotheses. Thus, type II error rates depend upon parameters of the 
distribution of the test statistic under the alternative hypothesis. The power of the test 
for a specified value of the hypothesized parameter is the probability the test would 
reject Ho when Ho is false, or 1 — p. 

In addition, whatever the test procedure is, rejecting Ho when Ho is not true is 
usually something we want a test to achieve with as high a probability as possible. 
Therefore we want the power of the test, for a given value of non-centrality parameter, 
to be high. Power algorithms are used to assist in the "best" selection of parameters 
for a statistical design. There do not exist simple closed form equations for computing 
the required probabilities from the common noncentral sampling distributions. 
Consequently, one must use numerical approximations to determine power. In this 
thesis we use various approximation methods to compute the cumulative probabilities 
for the noncentral F, noncentral x arid noncentral t distributions. 
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II. NONCENTRAL CHI-SQUARE, T,F DISTRIBUTIONS 



A. NONCENTRAL CHI-SQUARE DISTRIBUTION 

If \V1, \V2, .... Wn arc independently distributed as ,1), i= 1,2, ...,n. then 
Wi"^ has a distribution known as a noncentral X“ (n, X ) where n is the number of 
degrees of freedom and ^ M j ^ is the noncentrality parameter. When n j = H 2 
= ... = I* n ~ L = 0, and the noncentral (n,0) reduces to the usual central 

X'^ (n) with n degrees of freedom. The cumulative distribution function of x" (n,X ) is, 






- e 



i-A ^ V 



1 



}7i, j! p(i.j) 






(2.1) 

X X o 



while F(x;n, k) = 0 for x < 0 

It is possible to express F(x;n,X) for x>0, in an easily remembered form as a 
weighted sum of central x CDF's with weights equal to the probabilities of a Poisson 
distribution with expected value X/2. That is, 



p(Xl;r\,A) 

)=0 






( 2 . 2 ) 
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Thus a X (n, X) variable can be regarded as a mixture of central x variables. This 

2 

interpretation is often useful in deriving the distribution of functions of noncentral x 
random variables. 

The probability density function can, similarly, be expressed as a mixture of 
2 

central X probability density functions 




rtiu) 



where 




(2.3) 
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The mean and variance of the distribution are 

E(X)=n + 2X and Var(X)= 2n + 8X (Ref. l:p. 130]. 

B. NONCENTRAL F DISTRIBUTION 

IfYj and Y 2 are independent and Yj is x j ^ (Uj ,X ) and Y 2 is x ^ ( 1^2 ) then 



V = (Y,/n^)/(Y2/n2) 



(2.4) 



is distributed as F'(n^ ^2 »X ), the noncentral F distribution with n^ and n^ degrees of 
freedom and noncentrality parameter X . Its density function is given by 



CP 






where 



c = 



P (-^+Kj ( Tii + n, 
2*j 2U- >1, 



(2.5) 



when X =0 this reduces to the density function of the central F(npn 2 ) distribution. 
The mean and variance of the distribution are 






and 



Var(V) 



= r 



(i1ri)(i1i-+) 4- 



(for n, > 4) 



When X =0, these reduce to the mean and variance of the central F(nj,n 2 ) 
distribution. Derivation of (2.5) is shown in [Ref l:p. 189]. 

The cumulative distribution of V can be expressed in terms of an infinite series of 
multiples of incomplete beta functions, as follows; 
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is the incomplete beta function. 



where I yi, 

’''xf n,fo 




In this thesis, we use Paulson's Approximation from Severo and Zelen [Ref 2j. 
Their approximation is as follows: 



Pr(V < f 0 ) - (x) 



(2.7) 



where 

x = 



t 





'-1 


fi- 




.+A() ] 






•mV 1 




i- 



and 



0 is the standard Normal CDF. 

The following table shows some values obtained using this approximation, 
together with exact values of Pr(V < fo ). It can be seen in TABLE 1 that this 
approximation gives about 2 decimal place accuracy. Thus, this approximation 
provides adequate accuracy for computing the power of F-tests [Ref 3:p. 202]. 

The program we used in the computation of noncentral F CDF values is shown 
in Appendix D. 



C. NONCENTRAL T DISTRIBUTION 

The ratio 



T = (U + >,)/ V(Y/n) (2.8) 

where U and Y are independent random variables distributed as standard normal 
(N(0,1)) and x ^ with n degrees of freedom ,respectively, is said to have a noncentral t 
distribution with n degrees of freedom and noncentrality parameter X. Sometimes 
(or even .SX^), rather than X, is termed the noncentrality parameter. If X is equal to 
zero, the distribution is central t with n degrees of freedom. 

The cumulative distribution of T is [Ref l;p. 201] 
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Pr(T < t) =• 






P(^) 







(2.9) 



TABLE 1 

CDF OF A NONCENTRAL F DISTRIBUTION 



nl 


n2 


X 


fo 


approx 


exact 


3 


10 


4 


3.708 


0.7499 


0.745 






4 


6.552 


0.9192 


0.918 






16 


3.708 


0.2023 


0.206 






16 


6.552 


0.5190 


0.517 


3 


20 


4 


3.098 


0.7067 


0.700 






4 


4.938 


0.8894 


0.887 






16 


3.098 


0.1186 


0.126 






16 


4.938 


0.3488 


0.347 


5 


10 


6 


3.326 


0.7337 


0.731 






6 


5.636 


0.9143 


0.914 






24 


3.326 


0.1553 


0.158 






24 


5.636 


0.4629 


0.461 


5 


20 


6 


2.711 


0.6685 


0.664 






6 


4.103 


0.8715 


0.870 






24 


2.711 


0.0643 


0.069 






24 


4.103 


0.2437 


0.245 


8 


10 


9 


3.072 


0.7159 


0.714 






9 


5.057 


0.9087 


0.908 






36 


3.072 


0.1166 


0.119 






36 


5.057 


0.4088 


0.408 


8 


30 


9 


2.266 


0.581 


0.578 






9 


3.173 


0.8157 


0.813 






36 


2.266 


0.0146 


0.017 






36 


3.173 


0.0846 


0.088 



In this thesis, we use an approximation method suggested by Abamowitz and 
Stegun (Ref. 4:p. 949]. 

Pr(T < t) ~ a> (x) (2.10) 

where x = - and O is the standard Normal CDF. 
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The following table shows some values obtained using this approximation, 
together with exact values of Pr(T ^ to ). It can be seen in Table 2 that this 
approximation gives about 2 decimal place accuracy. Thus, this approximation 
provides adequate accuracy for computing the power of t-tests. 

The program we used in the computation of noncentral t CDF values is shown 
Appendix E. 



TABLE 2 

CDF OF A NONCENTRAL T DISTRIBUTION 



n 


X 


to 


approx 


exact 


4 


0.112 


4 


0.017 


0.01 




7.944 


4 


0.993 


0.99 




2.24 


10 


0.026 


0.01 




18.494 


10 


0.993 


0.99 


24 


5.786 


9.797 


0.012 


0.01 




13.981 


9.797 


0.991 


0.99 




16.118 


24.494 


0.014 


0.01 




33.078 


24.494 


0.992 


0.99 
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III. POWER OF THE F-TEST 



A. INTRODUCTION 

Under appropriate conditions, the best test for testing equality of several means 
is the analysis of variance test. Analysis of variance has a wide application. It is one 
of the most useful techniques in the field of statistical inference. As in any 
hypothesis-testing situation, the power of the F test is of interest to the experimenter. 

In this chapter we will discuss the power of F tests and provide an example. To 
give an overall evaluation of the power of F tests in the analysis of variance, we may 
use power curves. An important use of the power curve is to guide the experimenter in 
selecting fhe sample size (number of replicates) so that the design will be sufficiently 
sensitive to important potential difierences in the treatments. We will consider power 
curves for one-way and multi-way analyses of variance (AN'OVA's). 

B. THE ONE-WAY CLASSIFICATION ANALYSIS OF VARIANCE 

Suppose we have k different levels of a single factor that we wish to compare. 
The different levels of the factor are often called treatments. The observed responses 
from each of k treatments is a random sample on a random variable. The data would 
appear as in Figure 3.1. 



Treatment 



Observation 



1 


Yll 


Yi2 


■ 


2 


'^21 


Y 22 




• 


. 


. 


• 


k 


Ykl 


CNJ 

>- 


'*'kn 



Figure 3.1 Typical Data for One-way Classification ANOVA. 
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We will find it useful to describe the observations by the linear statistical model: 



Y- = H + x- + e- (3.1) 

where Yjj is the (ij)^^ observation, ^ is a parameter common to all treatments called 
the overall mean, tj is a parameter peculiar to the i^^ treatment called the i^^ treatment 
effect, and Cjj is a random error component, assumed to be I ID X (0, ). 

The indices used are: 

• i = the treatments, i = l,2,...,k 

• j = the replication per treatment, j = 1,2,. ..,n 

The objective in the AX'OVA is to test appropriate hypotheses about the 
treatment effects. The variance a~ is assumed constant for all levels of the factor. 
This model is called the one-way classification analysis of variance because only one 
factor is investigated. 

We are interested in testing the equality of the k treatment effects, so the 
appropriate hypotheses are 

flo • ~ ^2 ~ ■■■ ~ ^k ~ ^ 

Ha : Tj # Tj for some i,j . 

That is, if the null hypothesis is true, then each observation is made up of the mean ^ 
+ t plus a realization of the random error £”. 

The results of the ANOVA procedure is summarized in Table 3. 

We may use the expected values of the mean squares to verify that Fo (in Table 
3) is an appropriate test statistic for Ho. From the expected mean square we see that, 

9 

in general, MSE is an unbiased estimator of cr . Also, under the null hypothesis, MSt 

is an unbiased estimator of a . However, if the null hypothesis is false, then the 

expected value of MSt is +(n ^ (Tj)^)/(k-l). The expected value of the mean 
2 

square error is 

Therefore, under the alternate hypothesis the expected value of the numerator of 
the test statistic (Fo) is greater than the expected value of the denominator and we 
would reject Ho with values of the test statistic which are too large. That is, we would 
reject Ho if 

^o ^ ^a-(k-l)’(N-k) 
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TABLE 3 

THE ANALYSIS OF VARIANCE FOR A ONE-WAY CLASSIFICATION 



Source of 
Variation 


Sum of 
Squares 


Degrees of 
Freedom 


Mean 

Square 


Fo 


Between 

Treatments 


sst 


k-1 


MSt 


MSt/MSE 


Error(within 

treatments) 


SSE 


k(n-l) 


MSE 




Total 


SST 


kn-1 







where the numerator and denominator df are k-I and k(n-I) = N-k, and a is the type 1 
error rate. 

The power of the test is: 



1-p = P(Fq > I ^^Ise) 



(3.2) 



To evaluate the P in Equation 3.2 we need to know the distribution of the test statistic 
Fq if the null hypothesis is false. It can be shown that, if Ho is false, the statistic F^ 
has the noncentral F distribution with k-I and k(n-l) degrees of freedom and 
noncentrality parameter X, given by 

n I ( ^ i 

X = 



The noncentrality can be interpreted as the squared standardized distance 

between the origin and (tj,t 2 , ... The ratio ^ (t- is called the squared 

7 2 

standardized distance. If only an estimate of a is available, one may replace with 
the estimate [Ref. 5:p. 34). 



C. THE MULTI-WAY CLASSIFICATION ANALYSIS OF VARIANCE 

Many experiments require a study of the effects of two or more factors. It can 
be shown that, under certain conditions, factorial arrangements are the most efficient 
designs for this type of analysis. 
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One of the simplest factorial experiments involves only two factors or sets of 
treatments say factor A and factor B. Suppose there are a levels of factors A and b 
levels of factor B, and these are arranged in a 2-way factorial design; that is, each 
replication of the experiment contains all ab treatment combinations. Assume there 
are n replications of the experiment, and let represent the observation taken under 
the i''^ level of factor A and the level of factor B in the replication. 

The data can be summarized as shown in Figure 3.2. The order in which the abn 
observations are taken is selected at random. 



Factor B 





1 1 


1 1 


2 1 . . 


. 1 b 






1 1 


^iH’^112 1 


Y 121 .Y 122 1 


1 ^lbl»'^lb2 






1 1 1,. 


Yi2n 1 


’•••’ ^12n 1 

1 _ 


•• 1 >•••> Yibn 

_ l_ 






1 1 


Y 211 .Y 212 1 


^221 ’^222 1 


1 ^2bl’'^2b2 




Factor 


1 2 1 


>•••> ^21n 1 


, Y22n 1 
1 


• • 1 >• • • > Y2bn 

_ 1 




A 




... 1 


1 

1 


1 

1 






1 1 


Yall.Yal2 1 


Ya21>Ya22 1 


1 Yabl’^ab2 






1 3 1 


‘ > ^a2n 1 


«•••« Ya2n 1 • 

1 


• • 1 >• • • » Ygbn 

1 





Figure 3.2 Typical Data notation for a Two-way Classification.. 

The observations may be described by the linear model; 

Yijk = H + ti + pj+ Tpi- -HEi-k, (3.3) 

where 
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• [I = the overall mean efiect; 

• t j = the true efiect of i level of factor A, i = 1,2,. ..,a; 

• p j = the true efiect of j level of factor B, j = 1,2, ...,b; 

• r P •; = the effect of the interaction between t j and P ; ; and 

• ^ ijk “ ^ random error component, assumed to be I ID N (0,(T ). 

Both factors are assumed to be fixed. It is usually assumed that the treatment effects 

are defined as deviations from the overall mean, so ^ t j = 0 and X! P j ^ 

Similarly, the interaction effects are fixed and usally defined so that ^ t P ) jj = 0, 

where the summation is over either i or j. Since there are n replicates of the 

experiment, there are a total of abn observations. 

We are interested in testing various hypothesis about the parameters in equation 
3.3 . An appropriate hypothesis testing procedure would again be analysis of variance. 
.Vlore specifically, as we are considering two controllable sources of variation (A and 
B), the procedure is called the two-way classification analysis variance. 

In order to test the hypothesis Ho ; T j = 0 for i= 1,2,. ..,a (no row factor 

effects). Ho ; p j = 0 for j= l,2,...,b (no column factor effects), and Ho ; ( T p ) ” = 0 

(no interaction effects), we can express the total sum of square as; 

SST = SSA + SSB + SSAB + SSE. (3.4) 

Here, SSA is a sum of squares due to "rows" or factor A, SSB is a sum of squares due 
to "columns" or factor B, SSAB is a sum of squares due to the interaction between A 
and B, and SSE is a sum of squares due to error. The degrees of freedom associated 
with each sum of squares are shown in Figure 3.3. 

If we assume £ jjj^ are IID N(0,<y^ ) and apply Cochran's theorem (Theorem 3-1) 
under the null hypothesis of no effects, each sum of squares on the right-hand side of 
Table 4 when divided by c is distributed as x with the indicated number of degrees 
of freedom, and these statistics are independent. 

Theorem 3-1 (COCHRAN). Let Z j be HD N(0,1) for i= l,2,...,v and suppose 
j ^ = Qj + Q2 + Q3 +... + Q5 where s ^ v, and the quadratic form Qj has 

Vj degrees of freedom (i= l,2,...,s). Then Qj Q2 ,•••. Qj are independent 

chi-square random variable with Vj V2 degrees of freedom, respectively, if 

and only if 

V = Vj + V2 + ...+ Vg. 
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Effect 


Degrees of freedom 


A 


a-i 


B 


b-1 


AB interaction 


(a-D(b-l) 


Error 


ab(n-l) 


Total 


abn-1 



Figure 3.3 Table of degrees of freedom with sums of squares.. 

Assuming that factors A and B are fixed, the expected values of the 
mean squares are: 

E(MSA) = (T^ +(bnXTi^)/(a-l); (3.5) 

E(MSB) = ^ p . 2 y(b.i) . (3 6) 



E(MSAB) = + (n X I (t P ) ij ^ )/ (a-l)(b-l) ; and (3.7) 

E(MSE) = (P- . (3.8) 

Therefore, to test the hypothesis Ho : r j = 0 ; i = 1,2,. ..,a (no row factor effects), and 
Ho ; p j = 0 ; j = 1,2,.. .,b ( no column factor effects), and Ho : ( T P ) jj = 0 (no 
interaction effects), we would divide the corresponding mean square by the mean 
square error. Under the null hypothesis of no effect, this ratio will follow an F 
distribution with appropriate numerator degrees of freedom and ab(n-l) denominator 
degrees of freedom, and the critical region will be located in the upper tail. The test 
procedure is usually summarized in an analysis of variance table, such as shown in 
Table 4. 
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TABLE 4 

THE ANALYSIS OF VARIANCE FOR A TWO-WAY CLASSIFICATION 



Source of 
Variation 


Sum of 
Squares 


Degrees of 
Freedom 


Mean 

Square 


^0 


A treatments 


SSA 


a-1 


MSA 


MSA/MSE 


B treatments 


SSB 


b-l 


MSB 


MSB/MSE 


Interaction 


SSAB 


(a-D(b-l) 


MSAB 


MSAB/MSE 


Error 


SSE 


ab(n-l) 


MSE 




Total 


SST 


abn-1 







To the compute the power of tests in two-way ANOVAs. the procedure is the 
same as in the one-way case. Relationships among X, the numerator degrees of 
freedom and the denominator degrees of freedom are shown in Table 5 . 

In a similar way one can expand to the general multi-way ANOVA procedure 
[Ref 5:p. 124]. 

D. THE ALGORITHMS AND FLOWCHART 

A program for computing a power table and power curve in general ANOVA's 
(fixed model) is shown in Appendix F. The power of one-way ANOVA's are solved 
with these programs, using the following the sequence: 

1) Determine what variable to include, such as sample size vs power, number of 
treatments vs power, a-level vs power, or noncentrality vs power. 

2) Determine the values of the maximum, minimum, and increment of the variable 
chosen. 

3) Given input data, compute the critical value, using the inverse of the central 
F-distribution (Appendix A). 

4) Compute the power value in each case; that is, the CDF of the noncentral 
F-distribution (Chapter II). 

5) Print the power table or the power curve. 
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TABLE 5 

NONCENTRALITY PARAMETERS FOR POWER IN A TWO-WAY ANOVA 



Factor 


Numerator 

DOF 


Denomi nator 
DOF 


A 




a-1 


ab(n-l) 


B 




b-1 


ab(n-l) 


AB 


nZZ^Pjj^ 


(a-D(b-l) 


ab(n-l) 



Multi-way ANOVAs have the same algorithms but multi-way ANOVA may 
include tests of interaction effects. The program considers only up to three-way 
interaction effects. We can explain how to compute degrees of freedom of error term 
in m-way ANOVA's, assuming the balanced case, as follows: Total degrees of freedom 
is DOF(Total) = n ^ k(i) -1, where k(i) is the number of levels of the i*’*^ factor [Ref 9]. 

(1) If the model has only main effects without any interaction effects, 

DOFl(Error) = DOF(Total) - ^ (k(i)-l)- 

(2) If the model has several factors and only 2- way interaction effects, 

DOF2(Error)= DOFl(Error) - ^ S (k(i)-l)(k(j)-0- 

(3) IF the model has several factors and only 3-way interaction effects, 

DOF3(Error)=DOFl(Error) - S I Z (k(i)-l)(k(j)-l)(k(k)-l). 

(4) If the model has several factors and 2- way and 3- way interaction 

effects, 

DOF4(Error) = DOF3(Error) - Z Z (k(i)-l)(k(j)-l). 

If the model has more than 3-way interaction effects, then the user must modify 
the degrees of freedom of the error terms accordingly. A flowchart is shown in Figure 
3.4. 
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I'igure 3.4 System flowchart for F-test power. 
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E. EXAMPLES OF F-TESTS 



1. Example of One-way ANOVA test (number of replicates vs power) 

a. Scenario 

A manufacturer suspects that the batches of raw material furnished by his 
supplier differ significantly in calcium content. There are a large number of batches 
currently in the warehouse. Five of those are randomly selected for study. A chemist 
wants to know the appropriate sample size per each batch, in order to test 
Ho • ~ '** ~ ^k 

against 



Ha : t| 9^ Tj for some i,j 



The chemist would like to know how many replicates to run if it is important to reject 
Ho with probability at least 0.9 when the standardized distance square 
is 2 and a = .05. Thus he would like to know what the power of the F-test is for a 
range of possible replicates. He decides to check replicates from 2 to 12 in increments 
of 1. 

b. Inputs 

1) Select one-way ANOVA test (the number of replicates vs power). 

2) The number of treatment = 5. 

3) a -level = .05 

4) Standardized distance square = 2. 

5) M.^ximum replicates is 12, minimum replicates is 2, increment is 1. 

c. Start program 

Do you want to analyze one way ANOVA(y/n)? 

y 

Do you want to plot n (# of observation per treatment vs power)(y/n)? 

y 

Number of k( # of treatment)? 

7 

5 

Alpha- level? 

7 

.05 



Standardized squared distance value ? 

7 
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2. 

Maximum n value (the maximum value on X-axis)? 

? 

12 

Minimum n value (the minimum value on X-axis)? 
(condition :N must be more than 2) 

7 

2 

Increment n value ? 

7 



d. Output 

The output is presented in Table 6 

TABLE 6 

OUTPUT OF THE ONE-WAY ANOVA EXAMPLE (REPLICATES VS POWER) 



DOFl 


DOF2 


a 


F-INVERSE 


X 


POWER 


4 


5 


0.05 


5.29087 


4.00 


0.15373 


4 


10 


0.05 


3.52496 


6.00 


0.30024 


4 


15 


0.05 


3.10397 


8.00 


0.44935 


4 


20 


0.05 


2.91676 


10.00 


0.58677 


4 


25 


0.05 


2.81108 


12.00 


0.70332 


4 


30 


0.05 


2.74323 


14.00 


0.79562 


4 


35 


0.05 


2.69600 


16.00 


0.86449 


4 


40 


0.05 


2.66122 


18.00 


0.91328 


4 


45 


0.05 


2.63455 


20.00 


0.94631 


4 


50 


0.05 


2.61345 


22.00 


0.96776 


4 


55 


0.05 


2.59634 


24.00 


0.98118 



# of replicate POWER 



NN = 


2. 


POWER = 


0.15373 


NN = 


3. 


POWER = 


0.30024 


NN = 


4. 


POWER = 


0.44935 


NN = 


5. 


POWER = 


0.58677 


NN = 


6. 


POWER = 


0.70332 
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NN = 


7. 


POWER = 


0.79562 


NN = 


8. 


POWER = 


0.86449 


NN = 


9. 


POWER = 


0.91328 


NN = 


10. 


POWER = 


0.94631 


NN = 


11. 


POWER = 


0.96776 


NN = 


12. 


POWER = 


0.98118 



0. 98118 +-+ + + + + + — * *-+ 

POWER I * I 



0. 56746 



+ 






0. 15373 


+ 








1 

1 

1 

+ 

1 

1 

1 

1 

1 

1 

1 

1 

1 

+ 

1 

1 

1 

1 

1 

1 

1 

1 

1 

+ 

1 

+ 




2. 000 






7.000 


# OF REPLICATES 12.00 



X-SCALE: 0.1 25E + 00 UNITS 

Y-SCALE: "\"= 0.138E-01 UNITS 
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2. Example of one-way ANOVA test (noncentrality vs power) 

a. Scenario 

Five brands of batteries are under study. It is suspected that the life (in 
weeks) of the five brands is different. Five batteries of each brand are tested. A 
manufacturer wants to know the power as a function of the standardized squared 
distances. That is, would like to know what the power of the F-test is for a range of 
possible standardized distances square. He decides to check standardized squared 
distances from 1 to 5 in .2 increment, where n = 5 and a = .05. 

b. Inputs 

1) Select one-way ANOVA test (the standardized distance square vs power). 

2) Number of treatment = 5. 

3) Number of replicates = 5. 

4) a -level = .05. 

5) Maximum standardized distance square is 5 and minimum standardized distance 
square is 1 and the increment is .2. 

c. Start program 
Start program: 

Do you want to analyze one-way ANOVA(y/n)? 

y 

Do you want to plot n(# of observation per treatment) vs power(y/n)? 
n 

Do you want to plot alpha-level vs power(y/n)? 
n 

Do you want to plot noncentrality vs power(y/n)? 

y 

# of observations per treatment (n= ) ? 

7 

5 

# of treatment (k= )? 

7 

5 

Alpha-level ? 

7 

.05 

Maximum standardized squared distance value range ? 
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7 

5 

Minimum standardized squared distance value range ? 

7 

1 . 

Increment standardized squared distance value range ? 



.2 

d. Output 

The screen output ( Table 7 ) is as follows: 



TABLE 7 

OUTPUT OF THE ONE-WAY ANOVA EXAMPLE (NONCENTRALITY VS 

POWER) 



DOFl 


DOF2 


a 


F-INVERSE 




POWER 


4 


20 


0.05 


2.91676 


5.00 


0.30382 


4 


20 


0.05 


2.91676 


6.00 


0.36314 


4 


20 


0.05 


2.91676 


7.00 


0.42203 


4 


20 


0.05 


2.91676 


8.00 


0.47946 


4 


20 


0.05 


2.91676 


9.00 


0.53461 


4 


20 


0.05 


2.91676 


10.00 


0.58677 


4 


20 


0.05 


2.91676 


11.00 


0.63550 


4 


20 


0.05 


2.91676 


12.00 


0.68050 


4 


20 


0.05 


2.91676 


13.00 


0.72163 


4 


20 


0.05 


2.91676 


14.00 


0.75885 


4 


20 


0.05 


2.91676 


15.00 


0.79223 


4 


20 


0.05 


2.91676 


16.00 


0.82192 


4 


20 


0.05 


2.91676 


17.00 


0.84812 


4 


20 


0.05 


2.91676 


18.00 


0.87108 


4 


20 


0.05 


2.91676 


19.00 


0.89107 


4 


20 


0.05 


2.91676 


20.00 


0.90835 


4 


20 


0.05 


2.91676 


21.00 


0.92321 


4 


20 


0.05 


2.91676 


22.00 


0.93591 


4 


20 


0.05 


2.91676 


23.00 


0.94672 


4 


20 


0.05 


2.91676 


24.00 


0.95586 


4 


20 


0.05 


2.91676 


25.00 


0.96356 



,.-i\ 



29 



POWER 



NONCENTRALITY 



NN= 5.0 


POWER = 


0.30382 


NN= 6.0 


POWER = 


0.36314 


NN= 7.0 


POWER = 


0.42203 


NN = 8.0 


POWER = 


0.47946 


NN= 9.0 


POWER = 


0.53461 


NN = 10.0 


POWER = 


0.58677 


NN= 11.0 


POWER = 


0.63550 


II 

K) 

O 


POWER = 


0.68050 


NN= 13.0 


POWER = 


0.72163 


NN= 14.0 


POWER = 


0.75885 


NN= 15.0 


POWER = 


0.79223 


N’N= 16.0 


POWER = 


0.82192 


NN= 17.0 


POWER = 


0.84812 


NN= 18.0 


POWER = 


0.87108 


NN= 19.0 


POWER = 


0.89107 


NN= 20.0 


POWER = 


0.90835 


NN= 21.0 


POWER = 


0.92321 


NN= 22.0 


POWER = 


0.93591 


NN= 23.0 


POWER = 


0.94672 


NN= 24.0 


POWER = 


0.95586 


NN= 25.0 


POWER = 


0.96356 
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0.96356 

POWER 



+-+ 



+ + + + + * — 



* 



* 



0.63369 



+ 



* 



0.30382 


+-* +- 

5.000 


15.00 


NONCENTRALITY 




X-SCALE: ' 


•-"= 0.250E + 00 UNITS 






Y-SCALE: " 


r= O.IIOE-OI UNITS 
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IV. POWER OF THE T-TEST 



A. INTRODUCTION 

In testing the hypothesis that the mean of a normal distribution is equal to Pq 
when the standard deviation is unknown, the t statistic may be used. This test-statistic 
is a function of the sample mean, the sample standard deviation, and the sample size. 

The t-test has a very wide application in research areas. For example, we may 
wish to know whether product A is better than product B; or whether the outputs of 
machines C and D form a homogeneous mass of product; and so forth. To answer 
such questions we can employ the t-test which is one of the most useful techniques in 
the field of statistical inference. 

In such hypothesis-testing situations the power of t-tests is usually of interest to 
the experimenter. We may use the power curve to evaluate the power of the t-test. 
The required sample size might be determined by referring to a set of power curves. 
We will cf^nsider the power curves for the one-sample and two-sample cases. [Ref 10] 

B. ONE-SAMPLE T-TEST 

We assume that a random sample, XI, X2, X3, ..., Xn of size n is taken from a 
normal population with mean ^ and standard deviation <J. Suppose we wish to test the 
hypothesis that fi ^ against the alternative hypothesis that > fig . The 
procedure is to reject the hypothesis Ho : ^ ^ ftg if ((X- ^g)Vn)/S> t ^ , and to 
accept Ho if otherwise where X and S are the sample mean and sample standard 
deviation and where t is the 1 - percentage point of the Student t-distribution 
with n-1 degrees of freedom, i.e., 

P( Student t > t^j ) = a (4.1) 

With critical region ( th® hypothesis ft = fig is rejected with probability 

a when Ho is true. If the mean of the normal distribution is > Hg , then the 
power of the above test is 

P( noncentral t > t gj | 6 = (( ft j - ^g) Vn )/ <r) = 1- P (4.2) 
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where 6 is the noncentrality parameter for the noncentral t-distribution with n-1 
degrees of freedom [Ref 6:p .320]. 

C. TWO-SAMPLE T-TEST 

Suppose we are given two random samples from two normal populations XI, X2, 
..., Xn and Yl, Y2, .... Ym, where the X's and Y's are independent with mean values 
and common unknown variances = < 3 ^ = cr . Suppose we wish to 

test the null hypothesis Ho ; ^ jiy against the alternative hypothesis Ha that 

> Py. The procedure is to reject the hypothesis Ho if (X- Y) / V ( Sp ^ (I/n -I- 1/m)) 

> t ^ where t is a critical value of the Student t-distribution based on m+n— 2 

degrees of freedom. In this case take 6 = ( — ^y ) / cr V (I/n + 1/na) so the 

probability of rejecting the hypothesis when > jly is equal to 

P( noncentral t > t ^ | 8 ) = I- P (4.3) 

where 8 is the noncentrality parameter for the noncentral t-distribution with m-f-n-2 
degrees of freedom [Ref 6;p .321). 

D. TWO-SIDED T-TEST 

For tests where the alternative hypothesis specifies that )l # in the test of 
section B or y in the test of section C, where the direction can be either up or 

down, i.e, ^ > fig or < Hg and fi^^ > fly or fi^ < fly, we have what is known as a 
two-sided test. In this case the rejection procedure for the test of section B specifies 
that Ho is rejected if either ((X- fig)Vn)/S > t or ((X- fig)V ti)/S < - ^ 2 

and the power of the test is 

P( noncentral t > t gj /2 I 8 ) -H P( noncentral t < -t gj /2 I 8 ) = 1 - P (4.4) 
where 8 = ( fij -fig ) V n / cr [Ref 6:p .323]. 



E. THE ALGORITHMS AND FLOWCHART 

Listings of programs for generating the power table and curve for t-tests are 
given in the Appendix E. 

The following steps are used to compute the power: 

1) Determining whether one-sample or two-sample t-tests should be used. 
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2) After determining the test, give the input data, (n, a , fig , <r ) 

3) Determine whether the one-sided or two-sided test should be used. 

4) Compute the critical region (inverse t-values). (The computation method is 
provided in Appendix B) 

5) Compute the power corresponding to the chosen user option. (The 
computation method is discussed in Chapter II) 

6) Plot the power corresponding to each above option. 

A flowchart is shown in Figure 4.1a 



F. EXAMPLES OF T-TESTS 

1. Example of one sample one-sided t-test 

a. Scenario 

The supervisor in an electronic instrument plant is concerned with the 
repair time of an electronic equipment. Standards specify that the mean repair time 
must be at most 40 hours. For power computation, it is assumed that the variance of 
repairing time is 4.0. The appropriate hypothesis are 
Ho : fi > 40, 

Ha : fi < 40. 

The experimenter decides to use a sample of n = 25 observations and a = .05. He would 
like to know what the power of the t-test is for a range of possible population means. 
He decides to check means from 38 to 40 in .1 increments. 

b. Inputs 

INPUTS: step 1: one sample one-sided t-test 

step 2: fig = 40 
step 3: O’ = 2.0 
step 4: a = .05 
step 5: sample size = 25 

step 6: maximum fij =40., minimum fij =38, increment = .1 

c. Start program 
START PROGRAM: 



Do you want to test one sample t-test(y/n)? 

y 

Sample size ? 

7 



25 

Alpha-level ? 
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Figure 4.1a System flowchart for t-test power. 
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TWo sanrilfi T-test 



Input : sample size j 
alpha- level j 




Step ^ 




Figure 4.1b System flowchart for t-test power(cont'd. 
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Figure 4.1c System flowchart for t-test power(cont'd.) 
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? 



.05 

Mu-zero value ? 

? 

40 

Standard deviation value ? 

7 

2 

Do you want to use one-sided test(y/n)? 

y 

Do you want to test Ho : mu ^ muO vs Ha : mu > muO (y/n)? 
n 

Do you want to test Ho : mu ^ muO vs Ha : mu < muO (y/n)? 

y 

Maximum mul value (the maximum value on X-axis)? 

7 

40 

Minimum mul value (the minimum value on X-axis)? 

(condition : mu 1 must be less than muO ) 

7 

38 

Increment mul value ? 

7 

.1 

d. Output 

The screen output is as follows: 



39 



^*1 


power 




NN= 38.000 


POWER = 


0.99933 


NN= 38.100 


POWER = 


0.99849 


NN= 38.200 


POWER = 


0.99677 


NN= 38.300 


POWER = 


0.99346 


NN= 38.400 


POWER = 


0.98742 


NX= 38.500 


POWER = 


0.97706 


NN= 38.600 


POWER = 


0.96028 


XN= 38.700 


POWER = 


0.93465 


NN= 38.800 


POWER = 


0.89772 


XN= 38.900 


POWER = 


0.84755 


XX = 39.000 


POWER = 


0.78326 


XX = 39.100 


POWER = 


0.70559 


XX = 39.200 


POWER = 


0.61711 


XX = 39.300 


POWER = 


0.52203 


XX = 39.400 


POWER = 


0.42564 


XX = 39.500 


POWER = 


0.33355 


XX = 39.600 


POWER = 


0.25053 


XX = 39.700 


POWER = 


0.17996 


XX = 39.800 


POWER = 


0.12338 


XX = 39.900 


POWER = 


0.08062 
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0. 99933 
power 






+ 



+ - + 



- + + 

5«c i< 

* 

•k 



>+ 



k 



k 



k 



0. 53997 + 



0. 80616E-01 +-+ + + +--- 

38. 00 38. 95 

X-SCALE: 0.237E-0I UNITS 

Y-SCALE: 1"= 0.153E-01 UNITS 



k 



+ 



k 



k 



k 



k I 

k j 

+ + 

39.90 



2. Example of One sample two-sided t-test 
a. Scenario 

The supervisor in a rocket propellant plant is concerned with the burning 
rate of a rocket propellent. Standards specify that the mean burning rate must be 40 
inches per second. For power computation, it is assumed that o = 4.0. The 
appropriate hypotheses are 

Ho : = 40 

Ha : 40 
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The engineer decides to use a sample of n=25 observations and a =.05. He would like 
to know the power of the t-test is for a range of possible population means. He 
decides to check means from 38 to 42 in .1 increments. 

b. Inputs 



INPUTS: 



step I: one sample two-sided t-test 

step 2: = 40 

step 3: (T = 2.0 

step 4: a = .05 

step 5: sample size = 25 

step 6: maximum ^| =42, minimum =38, increment = .1 



c. Start program 
START PROGRAM: 

Do you want to test one sample t-test(y/n)? 



y 

Sample size ? 

7 



25 

Alpha-level ? 

7 

.05 

Mu-zero value ? 

7 



40 

Standard deviation value ? 

9 



2 

Do you want to use one-sided test(y/n)? 
n 

Do you want to use two-sided test(y/n)? 

y 

Maximum mul value (the maximum value on X-axis)? 

7 



42 

Minimum mul value (the minimum value on X-axis)? 
(condition : mul must be less than mu- zero) 



42 



9 



38 

Increment mul value ? 

9 

.1 

d. Output 

The screen output is as follows: 



POWER 



NN= 38.000 


POWER = 


0.99770 


NN= 38.100 


POWER = 


0.99525 


NN= 38.200 


POWER = 


0.99073 


NN= 38.300 


POWER = 


0.98279 


NN= 38.400 


POWER = 


0.96965 


NN= 38.500 


POWER = 


0.94910 


NN= 38.600 


POWER = 


0.91875 


NN= 38.700 


POWER = 


0.87639 


NN= 38.800 


POWER = 


0.82057 


NN= 38.900 


POWER = 


0.75109 


NN= 39.000 


POWER = 


0.66944 


NN= 39.100 


POWER = 


0.57882 


NN= 39.200 


POWER = 


0.48379 


NN= 39.300 


POWER = 


0.38975 


NN= 39.400 


POWER = 


0.30192 


NN= 39.500 


POWER = 


0.22459 


NN= 39.600 


POWER = 


0.16066 


NN= 39.700 


POWER = 


0.11146 


NN= 39.800 


POWER = 


0.07708 


NN= 39.900 


POWER = 


0.05692 


NN= 40.000 


POWER = 


0.05028 


NN= 40.100 


POWER = 


0.05687 


NN= 40.200 


POWER = 


0.07698 


NN= 40.300 


POWER = 


0.11131 


NN= 40.400 


POWER = 


0.16045 
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NX= 40.500 


POWER = 


0.22433 


NN= 40.600 


POWER = 


0.30162 


NN= 40.700 


POWER = 


0.38942 


NN= 40.800 


POWER = 


0.48344 


NN= 40.900 


POWER = 


0.57848 


NN= 41.000 


POWER = 


0.66912 


NN= 41.100 


POWER = 


0.75081 


NN= 41.200 


POWER = 


0.82034 


XN= 41.300 


POWER = 


0.87621 


XX = 41.400 


POWER = 


0.91862 


XX = 41.500 


POWER = 


0.94901 


XX = 41.600 


POWER = 


0.96959 


XX = 41.700 


POWER = 


0.98275 


XX = 41.800 


POWER = 


0.99070 


XX = 41.900 


POWER = 


0.99524 
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0.99770 

power 



0. 52398 



0. 50275E 



- + |. 

* 

* 



+ 



+ 



I i( ic i(ic I 

01 +-+ + + + + +-+ 

38.00 39.95 fij 41.90 

X-SCALE: 0.487E-01 UNITS 

Y-SCALE: 0.158E>01 UNITS 
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V. DESCRIPTION OF THE POWER PROGRAM 



A. PROCEDURE OVERVIEW 

The interactive program included in this thesis is written in FORTRAN 77. It 
was written for uce on an 1BM370 from IBM3278 terminal. It is an interactive 
program. The command 'NONCENT' is all that is required to start this program. The 
executive file name initializes the virtual machine environment and asks the user 
questions about the choice of program and compilation requirements. It then activates 
the selected program. Both ANOVA and TTEST are interactive programs which will 
be discussed in detail in later chapters. The output from the programs is presented on 
the terminal screen. 

B. INSTRUCTIONS FOR PROGRAM ACCESS 

To start the program, make sure you have loaded (NONCEN, ANOVA, and 
TTEST) on your disk. In C.MS (operating system mode), type 'NONCEN' and you will 
see: 

Please provide the FILENAME for your VS FORTRAN program. 

Now type the program name you want, for example 'TTEST', the response is; 

Do you need to compile your program ? (y/n) 

If you want to run, type 'Y' and your program will be loaded. 

A detailed description of the ANOVA (Analysis of Variance Test ) is given in 
Chapter III and of the t-test in Chapter IV. 

Following the screen output from the selected test, the user will be asked some 
questions about the output: 

Do you wish to BROWSE your output? (Y) 
n 

Print your output file? (Y) 

n 

Do you wish to XEDIT the program file? (Y/N) 
n 

Do you wish to run the program again? (Y) 
n 

Then return to CMS mode. 
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VI. SUMMARY AND CONCLUSIONS 



Power considerations are useful in the design and assessment of statistical tests. 
The calculation of power usually involves noncentral distributions for which tables of 
probability are not available. A search was made for algorithms approximating the 
noncentral t, F and x distributions. Algorithms giving sufficient accuracy and 
making efficient use of computer resources have been implemented in this thesis. 
Listings of the FORTRAN code for these implementations are included. 

An interactive program to compute and display power curves for several t-test 
and F-test situations has been developed. This program is user friendly and is 
described in this thesis. A listing of the program is provided. This program should be 
useful to researchers, experiment designers and statisticians. 
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APPENDIX A 

AN INVERSE CENTRAL F APPROXIMATION 



Let P be the central F complementary CDF, and let Q be the 
complementary standard normal CDF [Ref. 4: p. 947]. 

Define yp by Q(yp)=p, and Fp by P(Fp | Vj )=p. 

Then Fp ~ exp(2w) 
where w = wl-(w2)(w3) 

wl = (yp V (h+X))/h 

w2 = ( 1 /(V 2 -l)-l/(vi -D) 

w3 = (X+(5/6)+2/(3h)) 

h = 2/((l/(vi -1)+1/(V2 -1)) 

X=(yp 2 -3)/6 

The program we used in the computation of critical values of F-tests is shown 
below. 



* APPROXIMATION TO THE F- INVERSE DISTRIBUTION. * 

* CALCULATE F-INVERSE GIVEN DOFl ,DOF2 , ALPHA . * 

* Q( F-ALPHA |DOFl,DOF2 )= ALPHA. * 

* * 



REAL ALPHA ,NUMER , DENUM , LAM ,XP > YP , FINVER , A >B 
DATA NF1,NF2, ALPHA/2 >^,.01/ 

DATA CO ,C1 ,C2/2 .515517,0 . 802855,0.010528/ 

DATA D1,D2,D5/1.«2788, 0.189269, 0.001508/ 

B=NFl/2. 

A=NF2/2. 

C 

C CALCULATE NORMAL QUANTILE. 

C 

T=SQRT(ALOG(l./(ALPHA**2. ))) 

NUMER =C0 +C 1*T +C 2*T**2 

DENUM=1 . +D1*T+D25^T**2 . +D5*T**5 . 

XP=T-(NUMER/DENUM) 

C 

C CALCULATE THE INVERSE F-DISTRIBUTION. 

C 

LAM=(XPi«€2-5. )/6. 

H2= (l./(2.*B-l. ))-(!. /(2.*A-1. )) 

W4= (l./(2.*A-l. ))+(!. /(2.*B-1. )) 

Wl = (LAM*(5./6. )-(2./(5.*H)) ) 

H5=( XP«( H+LAM )** . 5 )/H 

H=W5-W2*W1 

FINVER=EXP(2*H) 

WRITE(5,100) NF1,NF2, ALPHA, FINVER 

100 FORMAK 'O' , ■ NF1= ' ,15, 'NF2= ' ,15, ' ALPHA=' ,F5.4, •FINVERSE = ' ,F6.5 ) 

STOP 

END 

The following table shows some values obtained using this approximation, 
together with exact values of critical values of the F distribution. It can be seen in 
Table 8 that this approximation provides adequate accuracy for the computation of the 
critical value of the F-tests. 
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TAB LI' 8 

APPROXIMATION AND HXACT CRITICAL VALULS OF I'-I USTS 



a 


nl 


n2 


approx 


exact 


.05 


4 


4 


6.581 


6.39 




4 


6 


4.601 


4.53 




4 


10 


3.525 


3.48 




10 


20 


2.35 


2.35 




20 


30 


1.932 


1.93 




6 


4 


6.359 


6.16 




10 


6 


4.112 


4.06 




20 


0 


3.944 


3.87 




30 


20 


2.041 


2.04 




30 


30 


1.841 


1.84 


.10 


4 


4 


4.150 


4.11 




4 


6 


3.218 


3.18 




4 


10 


2.648 


2.61 




10 


20 


1.939 


1.94 




20 


30 


1.668 


1.67 




6 


4 


4.046 


4.01 




10 


6 


2.952 


2.94 




20 


6 


2.864 


2.84 




30 


20 


1.739 


1.74 




30 


30 


1.607 


1.61 
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APPENDIX B 

THE INVERSE CENTRAL T APPROXIMATION 



We can use the asymptotic expansion for the Inverse t-CDF, as 
follows. Let Q be the complementary CDF of standard normal, and let P 
be the central t CDF [Ref. 4: p. 949]. 

If P(t i V ) = l-2p and Q(x ) = p 

then 

tp~ Xp + (g^ (Xp )/ v)+ (92 (Xp ) / )+(g 3 (Xp ) / )+( 94 (Xp ) / ) 

where (x)=. 25(x^ +x), 

92 (x)=(l/96)(5x^ +16x^ +3x), 

93 (x)=( l/384)(3x^ +19x^ +17x^ ~15x), and 

94 (x)=(l/92160)(79x^ +776x^ +1482x^ -1920x^ -945x). 

The program we used in the computation of the critical value of the 
t-test is shown below. 



* * 

* APPROXIMATION TO THE T- INVERSE DISTRIBUTION. * 

* CALCULATE T-INVERSE GIVEN DOF ALPHA. * 

* Q( T-ALPHA I DOF )= ALPHA. * 

* * 



C 

C set initial conditions. 

C 

REAL ALPHA ,NUMER,DENUM, LAM, XP,G1,G2,G3,X,TINVER 
INTEGER N 

DATA N, ALPHA/9,. 025/ 

DATA CO ,C1 ,C2/2 .515517,0 .802853 ,0 . 010328/ 

DATA D1,D2,D3/1. 432788, 0.189269,0.001308/ 

C 

C CALCULATE NORMAL QUANTILE. 

C 

ALPHA=.05 

T=SQRT( ALOG( l./( ALPHA**2. ) ) ) 

NUMER=CC+C1*T+C2*T**2 

DENUM=1 . +D1*T+D2*T**2 . +D3*T**3 . 

XP=T-(NUMER/DENUM) 

C 

C CALCULATE THE INVERSE T-DISTRIBUTION . 

C 

X=XP 

G1=.25*(X**3+X) 

G2=( 1./96. )*(5.*X**5+16.*X**3+3.*X) 

G3 = ( 1./384. )9t(3.*X**7+19.*X**3+17.*X**3-15.*X) 

G4=( 1./92160. )*( 79.*X**9+776.*X**7+1482.*X**5-1592.*X**3-945.*X) 

T P =X+ ( Gl/N } + ( G2/( N*N ) ) + 1 G3/ ( N**3 ) ) + ( G4/I N**4 ) ) 

TINVER=TP 

WRITE( 3,100) N, ALPHA, TINVER 

100 FORMATC '0* DOF=' ,I3,4X, ' ALPHA= ' ,F5 . 3 ,6X, *TINVERSE= ’ ,F9.5 ) 

STOP 

END 

The following table shows some values obtained using this 
approximation together with exact values of critical values of the t 
distribution. It can be seen in Table 9 that this approximation 
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provides adequate accuracy for computing the critical regions of the 
t-tests. 



TABLE 9 

APPROXIMATION AND EXACT CRITICAL VALUES OF T-TESTS 



DOF 




a ( Type 


I error ) 






. 05 


. 01 




exact 


approx 


exact 


approx 


2 


2. 919 


2. 867 


6. 964 


6. 475 


4 


2. 132 


2. 126 


3. 746 


3. 701 


6 


1. 943 


1. 942 


3. 142 


3. 130 


8 


1. 833 


1. 833 


2. 821 


2. 818 


10 


1. 812 


1. 823 


2. 764 


2. 761 


12 


1. 783 


1. 783 


2. 681 


2. 680 


14 


1. 761 


1. 761 


2. 624 


2. 624 


16 


1. 746 


1. 746 


2. 583 


2. 583 


18 


1. 734 


1. 734 


2. 552 


2. 552 


20 


1. 725 


1. 725 


2. 528 


2. 528 


25 


1. 708 


1. 708 


2. 485 


2. 485 


35 


1. 690 


1. 690 


2. 438 


2. 438 


45 


1. 690 


1. 690 


2. 412 


2. 413 


65 


1. 669 


1. 669 


2. 385 


2. 385 


85 


1. 663 


1. 663 


2. 371 


2. 371 
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APPENDIX C 

AN APPROXIMATION OF THE NORMAL CDF 



Let Z be the standard normal density function and let P be the 
corresponding CDF [Ref. 4: p. 939]. Then 

P(X) = J Z(t)dt ~ l-Z(x)(ajt+a2t^+a3^3)+ c (x) 

where t=l/(l+px) 

1 € (x)| < .000001 

p =. 33267, and a^ = . 4361836 

a2 = - . 1201676 
a3 = .9372980 



52 



APPENDIX D 

THE COMPUTATION OF NONCENTRAL F DISTRIBUTION 



***************^«t****** INFORMATION ********************************** 



* * 

* THE OBJECTIVE OF THIS PROGRAM IS TO CALCULATE NON-CENTRAL F * 

* DISTRIBUTION ( CDF , F-INVERSE ) . * 

* THE CALCULATION METHOD IS NORMAL APPROXIMATION. * 

* * 

********************** VARIABLE DEFINITION *************************** 
* * 

* FVAL : X-VALUE DIMENSION * 

* FCDF : NON-CENTRAL F CDF DIMENSION. * 

* FINVER : GIVEN X VALUE, FIND INVERSE VALUE. * 

* FLAM : NON-CENTRAL PARAMETER * 

* FX : TO FIND INVERSE, GIVEN X VALUE. * 

* FSTAR : X VALUE VARIABLE. * 

* NFl : DEGREE OF FREEDOM 1. * 

* NF2 : DEGREE OF FREEDOM 2. * 

* CDF ; NORMAL DISTRIBUTION FUNCTION. * 

* P ( Z <= X) * 

* CONST : 1 / SORT ( 2 * PI ) IN NORMAL DISTRIBUTION FUNCTION. * 

* B : THE VALUE OF *X' IN NORMAL DISTRIBUTION. * 

* * 



************************************************************************ 

C 

C SET INITIAL CONDITIONS. 

C 

REAL FVAL( 100 ) , FCDF( 100 ) , FINVER , FLAM , FMAX , FX , FINI AL 
DATA NFl ,NF2 , FLAM , FMAX/2 ,<^,75.0,90./ 

DATA Al, A2, A3, P/. <^361836, -.1201676,. 9372980, .33267/ 

DATA FSTAR/2S.128/ 

HRITE(3,2^) NFl, NF2, FLAM 
HRITE(3,21) 

WRITE(3,22) 

HRITE(3,23) 

CONST = 1.0 / SQRT (2.0 * 3.1^15927) 

C 

C CALCULATE THE NORMAL APPROXIMATION . 

C 

Cl=( NFl^FSTAR )/( NFl+FLAM ) 

C2=l-(2./(9.*NF2)) 

C3=(2.*(NF1+2.*FLAM))/(9.*(NF1+FLAM)**2. ) 

CCl=(Cl**(l./3. ))*C2-( 1-C3) 

CC2 = ( (C3+(2./( 9.*NF2))*(Cl)**(2./3. ))**.S) 

FX=CC1/CC2 
CDF = 0.0 
IF(FX .EQ.O) THEN 
CDF = 0.5 

ELSE IF(FX.GT.O) THEN 
T=1./(1.+P*FX) 

XSQU=( FX*FX)*.S 
ZPDF=CONST*EXP( -XSQU ) 

CDF=1-ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

ELSE 

FX=-FX 

T=1./(1.+P*FX) 

XSQU=(FX*FX)*.S 
ZPDF=CONST*EXP( -XSQU ) 

CDF=ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

END IF 

FVAU FSTAR )=FSTAR 
FCDF( FSTAR )=CDF 

WRITE (3,250) FSTAR ,, FCDF ( FSTAR ) 

250 FORMAT( *0 ‘ ,F7.3,8X,F8.5 ) 

21 FORMAT! *0' ,40( )) 

22 FORMAT! 'O' ,2X,*"X"-VALUE' ,8X,'F(Z <= X)') 

23 FORMAT! 'O', 40! '_' )) 

24 FORMAT! '0 ', 'D0F1= ' ,12, 3X, 'DOF2=' ,12, 3X, 'NON-CENTRAL PARAMETER^ 
&',F8.4) 

STOP 

END 
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APPENDIX E 

THE COMPUTATION OF NONCENTRAL T DISTRIBUTION 



XXX X'XXXXX INFORMATION ^fXXX'XX'XXXXXXXXXXXXX'XXXX'XXXXXXXX'XX'X 



* THE OBJECTIVE OF THIS PROGRAM IS TO CALCULATE NON-CENTRAL T * 

* DISTRIBUTIONS C CDFs ) 

* THE CALCULATION METHOD IS NORMAL APPROXIMATION. * 

* (REF. SEVERO & ZELEN'S (1960) ) * 

* 

****9(^(*9He*****9t****9t*9t VARIABLE DEFINITION J^XXXXXXXXXX XXXXXXXXXXX X X 

X X 

* TVAL : X-VALUE DIMENSION * 

* TCDF : NON-CENTRAL T CDF DIMENSION. 

* TLAM : NON-CENTRALITY PARAMETER * 

* TX : TO FIND INVERSE, GIVEN X VALUE. * 

* TSTAR ; X VALUE VARIABLE. * 

* NF : DEGREE OF FREEDOM * 

* CDF : NORMAL DISTRIBUTION FUNCTION. * 

* P ( Z <= X) ^ 

* CONST ; 1 / SQRT ( 2 * PI ) IN NORMAL DISTRIBUTION FUNCTION. * 

^ B : THE VALUE OF 'X' IN NORMAL DISTRIBUTION. * 

* ^ 



c 

c SET INITIAL CONDITIONS. 

REAL TVAL( 1DD),TCDF( 100 ) ,TINVER , TLAM, TMAX,TX,TINIAL 
DATA NF,TLAM/2<h,33.D78/ 

DATA A1 , A2 , A3 ,P/ . <h361836 ,- . 1201676 , . 9372980 , . 33267/ 
DATA TSTAR/24.<h949/ 

WRITE(3,24) NF,TLAM 
HRITE(3,21) 

WRITE(3,22) 

WRITE(3,23) 

CONST = 1.0 / SQRT (2.0 * 3.1415927) 

CALCULATE THE NORMAL APPROXIMATION . 

Cl=( 1 . -( 1 . /( 4 . *NF ) ) )*TSTAR-TLAM 
C2=SQRT( 1 . +( TSTAR^t*2/( 2 . *NF ) ) ) 

TX=C1/C2 

CDF = 0.0 
IF(TX .EQ.O) THEN 
CDF = 0.5 

ELSE IF(TX.GT.D) THEN 
T=1./(1.+P*TX) 

XSQU=(TX*TX)*.5 
ZPDF=CONST*EXP( -XSQU ) 

CDF=1-ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

ELSE 
TX=-TX 

T=1./(1.+P*TX) 

XSQU=(TX*TX)*.5 
ZPDF=CONST*EXP( -XSQU) 

CDF=ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

END IF 

TVALC TSTAR )=TSTAR 
TCDF( TSTAR )=CDF 

WRITE (3,250) TSTAR ,, TCDF ( TSTAR ) 

250 FORMATC '0 ' ,F7. 3,8X,F8.5 ) 

21 FORMATC '0‘,40( *_' )) 

22 FORMATC '0‘,2X,*"X”-VALUE' ,8X,*F(Z <= X)') 

23 FORMATC *0*,4D( •_' )) 

24 FORMATC 'D* , 'DOF= ' ,12, 3X, 'NON-CENTRAL PARAMETER= 

S* ,F8.4) 

STOP 

END 
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APPENDIX F 

THE POWER OF ANOVA PROGRAM LIST 



* POWER ALGORITHM ( ANALYSIS OF VARIANCE ) * 

* * 

* DIRECTED BY ; PROFESSOR DONALD .R. BARR * 

* WRITTEN BY : HUR, SEONG PIL JULY 1986 * 

* DEPARTMENT OF O.R, NAVAL POSTGRADUATE SCHOOL * 

* * 

* PROGRAM USES THE POWER PLOT TO TEST FOR ANALYSIS OF VARIANCE. * 

* (A) ONE-WAY ANOVA * 

* (1) SAMPLE SIZE VS POWER * 

* (2) NUMBER OF TREATMENT VS POWER * 

* (3) ALPHA-LEVEL VS POWER * 

* (^) NONCENTRALITY PARAMETER VS POWER * 

* (B) MULTI-WAY ANOVA * 

* (1) NUMBER OF REPLICATES VS POWER * 

* (3) ALPHA-LEVEL VS POWER * 

* (^) NONCENTRALITY PARAMETER VS POWER * 

* NOTE ; subroutine PLOTT is NONIMSL subroutine library. * 

* * 

REAL NN( 200 ) , POWER! 200 ) ,NON , FINV ,NONCEN , FPOW , FX , ALPHA ,CDF , CONST 
REAL MAXDEL ,MINDEL ,INCDEL , STDS, PICON ,MINSD ,MAXSD ,INCSD ,MAXALP 
REAL MINALP,INCRAL 

INTEGER MF( 20 ) ,MAXN, MINN, INCRN ,MAXK , MINK ,INCRK,JJ ,K,N 
CHARACTER*! ANS 



START ONE-WAY OR MULTI-WAY ANALYSIS OF VARIANCE. 



00 PRINT*, *DO YOU WANT TO USE ONE-WAY ANOVA (Y/N)?* 

READ! 5,7) ANS 
7 FORMAT (Al) 

IF! !ANS.EQ. 'N* ).OR. ! ANS.EQ. 'N* ) ) GO TO 1000 

PRINT*, 'DO YOU WANT TO PLOT N! n OF OBSERVATION PER TREATMENT) 

& VS POWER ! Y/N)?' 

READ! 5,7) ANS 

IF! ! ANS.EQ. 'Y' ). OR. !ANS.EQ. *Y' ) ) GO TO 10 

PRINT*, 'DO YOU WANT TO PLOT K ! NUMBER OF TREATMENT ) VS POWER 
& ! Y/N)?' 

READ!5,7) ANS 

IF! !ANS.EQ. *Y' ).OR.!ANS.EQ. *Y' )) GO TO 20 

PRINT*, 'DO YOU WANT TO PLOT ALPHA-LEVEL VS POWER !Y/N)?' 

READ! 5,7) ANS 

IF! ! ANS.EQ. 'Y' ). OR. !ANS.EQ. 'Y' )) GO TO 30 

PRINT*, 'DO YOU WANT TO PLOT NONCENTRALITY VALUE VS POWER !Y/N)?' 
READ!5,7) ANS 

IF! !ANS.EQ. 'Y' ).OR.!ANS.EQ. 'Y' )) GO TO ^0 
PRINT*, 'DO YOU WANT TO PLOT AGAIN !Y/N)?' 

READ!5>7) ANS 

IF(!ANS.EQ. 'Y' ).OR.!ANS.EQ. 'Y' )) GO TO 500 
GO TO 999 



ANALYZE NUMBER OF SAMPLE SIZE PER TREATMENT VS POWER IN ONE-WAY 
ANALYSIS OF VARIANCE. 



10 PRINT*, 'NUMBER OF K! n OF TREATMENT )?' 

READ !5,*) K 
PRINT*, 'ALPHA- LEVEL?' 

READ !5,*) ALPHA 

PRINT*,' STANDARDIZED SQUARED DISTANCE VALUE ?' 
READ !5,*) STDS 
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PRINT*,’ MAXIMUM N VALUE ( THE MAXIMUM VALUE ON X-AXIS )?’ 

READ (5,*) MAXN 

PRINT*, ‘MINIMUM N VALUE ( THE MINIMUM VALUE ON X-AXIS )? 

S (CONDITION :N MUST BE MORE THAN 2 )' 

READ (5,*) MINN 

PRINT*, ‘INCREMENT N VALUE ?‘ 

READ (5,*) INCRN 
HRITE(6,95) 

95 FORMAT! IX, ‘DOFl DOF2 ALPHA-LEVEL F-INVERSE NONCENTRAL POWER’) 
NF1=K-1 
JJ=1 

DO 110 N=MINN, MAXN, INCRN 
NF2=K*(N-1) 



COMPUTE CENTRAL F-INVERSE GIVEN DOFl,DOF2 AND ALPHA-LEVEL. 



FINV=FINVER ( NFl ,NF2, ALPHA ) 



COMPUTE NONCENTRALITY PARAMETER AS SAMPLE SIZE CHANGE. 



NONCEN=N*STDS 



COMPUTE POWER USING NONCENTRAL F-DISTRIBUTION (CDF) GIVEN DOFl, 
D0F2, CENTRAL F-INVERSE AND NONCENTRALITY VALUE. 



POW=FPO( ALPHA ,NF1 ,NF2 , FINV ,NONCEN ) 

WRITE (6,15) NFl ,NF2 , ALPHA , FINV,NONCEN ,POW 
15 FORMAT! ‘ 0 ’ ,I5,3X,I5 ,2X, 2F10 .5 , IX, F12 . 2 ,1X, FID .5 ) 
NN( JJ)=N 
POWER! JJ)=POW 
JJ=JJ+1 
110 CONTINUE 
JJ=JJ-1 

PRINT*, 'NN POWER’ 

DO 76 N=1,JJ 

WRITE! 6 ,75 ) NN! N ) , POWER! N ) 

75 FORMAT! ‘O’, ’NN= ’ ,F8.3,^X, ’POWER= ’ ,F10.5 ) 

76 CONTINUE 



PLOT NUMBER OF SAMPLE SIZE PER TREATMENT VS POWER IN ONE-WAY 
ANALYSIS OF VARIANCE USING LIBRARY SUBROUTINE. 



CALL PLOTT! NN, POWER, JJ,0) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN !Y/N)?‘ 
READ (5,7) ANS 

IF! (ANS.EQ. ‘Y’ ).OR.!ANS.EQ. ’Y’ )) GO TO 500 
GO T:^ 999 



ANALYZE NUMBER OF TREATMENT VS POWER IN ONE-WAY ANALYSIS OF VARIANCE 



20 PRINT*, OF OBSERVATIONS PER TREATMENT (N= ) ?’ 

READ !5,*) N 

PRINT*, ‘ALPHA LEVEL (REAL VALUE) ?’ 

READ !5,*) ALPHA 

PRINT*, ‘WHAT STANDARDIZED SQUARED DISTANCE VALUE ?’ 

READ !5,*) STDS 

PRINT*, ‘MAXIMUM K !^^ OF TREATMENTS) RANGE ?’ 

READ (5,*) MAXK 

PRINT*, 'MINIMUM K (# OF TREATMENTS) RANGE !K MUST BE MORE THAN 2)?’ 

READ C5,*) MINK 

PRINT*, 'INCREMENT K VALUE ?’ 

READ (5,*) INCRK 
WRITE(6,195) 

195 FORMAT! IX, ‘DOFl D0F2 ALPHA-LEVEL F-INVERSE NONCENTRAL POWER’) 
JJ=1 

DO 111 K=MINK, MAXK, INCRK 

NF1=K-1 

NF2=K*(N-1) 

FINV=FINVER ! NFl, NF2 , ALPHA ) 

NONCEN=N*STDS 

POW=FPO! ALPHA ,NF1 ,NF2 ,FINV,NONCEN ) 

WRITE! 6,115) NFl ,NF2 , ALPHA , FINV ,NONCEN ,POW 
115 FORMAT! 'O' ,I5,3X,I5,2X,2F10 . 5,1X,F12 . 2 , IX, FID .5 ) 

NN!JJ)=K 
POWER! JJ )=POW 
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JJ=JJ+1 
111 CONTINUE 

PRINT3^,'NN POWER' 

DO 176 N=1,JJ 

WRITE( 6,175) NN(N),POWER(N) 

175 FORMAK'O', *NN=' ,F8.3,^X, 'POWER=' ,F10.5 ) 

176 CONTINUE 

CALL PLOTT (NN, POWER, JJ,0) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN (Y/N)?' 

READ (5,7) ANS 

IFHANS.EQ.'Y' ).OR.(ANS.EQ.'Y' )) GO TO 500 
GO TO 999 

C 

C 

C ANALYZE ALPHA-LEVEL VS POWER IN ONE-WAY ANALYSIS OF VARIANCE. 

C 

C 

30 PRINT*,'# OF OBSERVATIONS PER TREATMENT (N= ) ?' 

READ (5,*) N 

PRINT*,'# OF TREATMENTS ( K= )?' 

READ (5,*) K 

PRINT*,' STANDARDIZED SQUARED DISTANCE VALUE ?' 

READ (5,*) STDS 

PRINT*, 'MAXIMUM ALPHA RANGE ?' 

READ (5,*) MAXALP 

PRINT*, 'MINIMUM ALPHA RANGE ?' 

READ (5,*) MINALP 

PRINT*, 'INCREMENT ALPHA VALUE ?' 

READ (5,*) INCRAL 
WRITE(6,295) 

295 FORMATCIX, 'DOFl D0F2 ALPHA-LEVEL F-INVERSE NONCENTRAL POWER') 
JJ=1 
NF1=K-1 
NF2=K*(N-1) 

DO 211 ALPHA=MINALP, MAXALP, INCRAL 
FINV=FINVER ( NFl ,NF2 , ALPHA ) 

NONCEN=N*STDS 

POW=FPO( ALPHA, NFl, NF2,FINV,N0NCEN) 

WRITE! 6,215) NFl ,NF2, ALPHA, FINV,NONCEN,POW 
215 FORMAT! ' 0 ' ,I5,3X,I5, 2X ,2F10 .5 ,1X, F12 . 2 ,1X ,F10 .5 ) 

NN( JJ)=ALPHA 
POWER! JJ)=POW 
JJ=JJ+1 
211 CONTINUE 
JJ=JJ-1 

PRINT*, 'NN POWER' 

DO 276 I =1,JJ 

WRITE ! 6 ,275 ) NN! JJ ), POWER! JJ ) 

275 FORMAT! 'O', 'NN= ' ,F8. 3 ,^X, ' POWER= ' , FIO .5 ) 

276 CONTINUE 

CALL PLOTT! NN, POWER, JJ,0) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN !Y/N)?' 

READ !5,7) ANS 

IF!!ANS.EQ. 'Y' ).OR.!ANS.EQ. 'Y' )) GO TO 500 
GO TO 999 

C 

C 

C ANALYZE NONCENTRALITY PARAMETER VS POWER IN ONE-WAY ANALYSIS OF 

C VARIANCE 

C 

C 

40 PRINT*,'# OF OBSERVATIONS PER TREATMENT ! N= ) ?' 

READ !5,*) N 

PRINT*,'# OF TREATMENTS! K=)?' 

READ !5,*) K 

PRINT*, 'ALPHA-LEVEL ?' 

READ !5,*) ALPHA 

PRINT*, 'MAXIMUM STANDARDIZED SQUARED DISTANCE VALUE RANGE ?' 

READ (5,*) MAXSD 

PRINT*, 'MINIMUM STANDARDIZED SQUARED DISTANCE VALUE RANGE ?' 

READ .!5,*) MINSD 

PRINT*, 'INCREMENT STANDARDIZED SQUARED DISTANCE VALUE RANGE ?' 

READ (5,*) INCSD 
WRITE! 6, 495) 

495 FORMAT! IX, 'DOFl D0F2 ALPHA-LEVEL F-INVERSE NONCENTRAL POWER') 
JJ=1 
NF1=K-1 
NF2=K*!N-1) 

DO 411 STDS=MINSD, MAXSD, INCSD 
FINV=FINVER ! NFl ,NF2 , ALPHA ) 

NONCEN=N*STDS 

POW=FPO! ALPHA, NFl, NF2,FINV,N0NCEN) 

WRITE! 6,415) NFl ,NF2 , ALPHA, FINV,NONCEN,POW 
415 FORMAT! 'O' ,15 ,3X,I5,2X, 2F10 .5 ,1X, F12 . 2 , IX, FIO . 5 ) 

NN! JJ)=NONCEN 
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POWER( JJ)=POW 
JJ=JJ+1 
411 CONTINUE 
JJ=JJ-1 

DO 476 I =1,JJ 

WRITE ( 6 ,475 ) NN( I ) ,POWER( I ) 

475 FORMAT(*0', ’NN= ' ,F8.3,4X, •POWER= ’ ,F10 .5 ) 

476 CONTINUE 

CALL PLOTKNN, POWER, JJ, 0 ) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN (Y/N)?' 

READ (5,7) ANS 

IF((ANS.EQ. 'Y* ).OR.(ANS.EQ. ’Y' ) ) GO TO 500 
GO TO 999 

C 

C 

C ANALYZE MULTI-WAY ANALYSIS OF VARIANCE. 

C 

c 

1000 PRINT*, 'HOW MANY WAYS ANOVA CLASSIFICATION?’ 

READ(5,*) M 
DO 600 1=1, M 

PRINT*, 'FOR', I, '-FACTOR LEVEL ,HOW MANY LEVELS ?' 

READ(5,*) MF(I) 

600 CONTINUE 

C 

c 

C ANALYZE MODEL WITH FACTOR-LEVEL IN MULTI-WAY ANOVA. 

C 

c 

PRINT*, 'IS THERE ANY INTERACTION ?(Y/N)‘ 

READ(5,7) ANS 

IF( (ANS.EQ. 'Y' ).OR.(ANS.EQ. 'Y' )) GO TO 3100 

ISUM=0 

ITOT=l 

DO 800 I=1,M 

ISUM=ISUM+(MF(I)-1) 

ITOT=ITOT*MF(I) 

800 CONTINUE 

IERDF=ITOT-ISUM-l 

C 

c 

C COMPUTE MINIMUM NUMBER OF REPLICATE PER TREATMENT. 

C 

C 

MINR=MINIR( ITOT ,ISUM ) 

PRINT*, 'AT LEAST MINIMUM NUMBER OF REPLICATE PER TREATMENT =',MINR 
PRINT*, 'WHICH FACTOR ANALYSIS (IF YOU WANT 3 RD FACTOR, PRESS 3)?' 
READ(5,*) I 
NF1=MF(I)-1 

PICON=FLOAT( ITOT )/FLOAT( MF( I ) ) 

GO TO 2000 

C 

c 

C ANALYZE MODEL WITH FACTOR-LEVEL AND 2-WAY INTERACTION IN 

C MULTI-WAY ANOVA. 

C 

c 

3100 PRINT*, 'ARE THERE ONLY 2-WAY INTERACTION AND FACTOR-LEVEL?( Y/N ) ' 
READ(5,7) ANS 

IF( (ANS.EQ. *N' ). OR. (ANS.EQ. *N' )) GO TO 2150 

ISUM1=0 

ISUM2=0 

ITOT=l 

DO 850 1=1, M 

ISUMl=ISUMl+( MF( I )-l ) 

ITOT=ITOT*MF(I) 

850 CONTINUE 

DO 860 I=1,M-1 
DO 870 J=I+1,M 

ISUM2=ISUM2+( MF( I )-l )*( MF( J )-l ) 

870 CONTINUE 
860 CONTINUE 

ISUM=ISUM1+ISUM2 
IERDF=ITOT-ISUM-l 
MINR=MINIR( ITOT ,ISUM ) 

PRINT*, 'AT LEAST MINIMUM NUMBER OF REPLICATE PER TREATMENT =',MINR 
PRINT*, 'DO YOU WANT TO USE ONLY FACTOR-LEVEL (Y/N)?' 

READ(5,7) ANS 

IF( (ANS.EQ. 'N' ).0R.( ANS.EQ. 'N' ) ) GO TO 880 

PRINT*, 'WHICH FACTOR ANALYSIS (IF YOU WANT 3 RD FACTOR, THEN 
aPRESS 3)?' 

READ(5,*) I 
NF1=MF(I)-1 

PICON=FLOAT( ITOT )/FLOAT( MF( I ) ) 

GO TO 2000 
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880 PRINT*, ’WHICH INTERACTION ANALYSIS?( IF YOU WANT 2 AND 3 WAY INTER 
FACTION, THEN PRESS 2 3 ).* 

READ(5,*) I,J 
NF1=(MF(I)-1)*(MF( J)-l ) 

PICON=FLOAT( ITOT )/FLOAT( MF( I )*MF( J ) ) 

GO TO 2000 

C 

C 

C ANALYZE MODEL WITH FACTOR-LEVEL AND 3-WAY INTERACTION IN 

C MULTI-WAY ANOVA. 

C 

C 

2150 PRINT*, 'ARE THERE ONLY 3-WAY INTERACTION AND FACTOR-LEVEL ?(Y/N)' 
READ(5,7) ANS 

IF( (ANS.EQ. 'N* ).OR.(ANS.EQ. 'N' )) GO TO 2200 

ISUM1=0 

ISUM3=0 

IT0T=1 

DO 871 1=1, M 

ISUMl=ISUMl + ( MF( I )-l ) 

ITOT=ITOT*MF(I) 

871 CONTINUE 

DO 872 1=1, M-2 

DO 873 J=I+1,M-1 
DO 874 K=J+1,M 

ISUM3=ISUM3+( MF( I )-l )*( MFC J )-l )*( MFC K )-l ) 

874 CONTINUE 
873 CONTINUE 

872 CONTINUE 
ISUM=ISUM1+ISUM3 
IERDF=IT0T-ISUM-1 
MINR=MINIRC ITOT , ISUM ) 

PRINT*, 'AT LEAST MINIMUM NUMBER OF REPLICATE PER TREATMENT =*,MINR 
PRINT*, 'DO YOU WANT TO ANALYSIS ONLY FACTOR-LEVEL ?(Y/N)* 

READC5,7) ANS 

IFC CANS. EQ. 'N' ). OR. CANS. EQ. 'N* )) GO TO 875 

PRINT*, ’WHICH FACTOR ANALYSIS (IF YOU WANT 3 RD FACTOR, PRESS 3)?* 

READC5,*) I 

NF1=MFCI)-1 

PICON=FLOATC ITOT )/FLOAT( MFC I ) ) 

GO TO 2000 

875 PRINT*, ’DO YOU WANT TO ANALYSIS 3-WAY INTERACTION ?CY/N)’ 

READ(5,7) ANS 

IFCCANS.EQ. ’N’ ).OR.CANS.EQ. ’N’ )) GO TO 500 

PRINT*, ’WHICH INTERACTION ANALYSIS?C IF YOU WANT 2 AND 3 AND 4 INT 
&ERACTION ,THEN PRESS 2 3 4 

READC5,*) I,J,K 

NF1 = CMFC I )-l )*CMF( J)-1)*CMFCK)-1) 

PICON=FLOATC ITOT )/FLOATC MFC I )*MFC J )*MFC K ) ) 

GO TO 2000 

C 

C 

C ANALYZE MODEL WITH FACTOR-LEVEL ,2-WAY INTERACTION IN AND 

C 3-WAY INTERACTION IN MULTI-WAY ANOVA. 

C 

C 

2200 PRINT*, ’ARE THERE 2-WAY INTERACTIONS AND 3-WAY INTERACTIONS AND 
SFACTER-LEVEL (Y/N)?’ 

READC5,7) ANS 

IFC CANS. EQ. ’N’ ). OR. CANS. EQ. *N’ )) GO TO 999 

ISUM1=0 

ISUM2=0 

ISUM3=0 

IT0T=1 

DO 950 I=1,M 

ISUM1=ISUM1+C MFC I )-l ) 

ITOT=ITOT*MFCI) 

950 CONTINUE 

DO 960 1=1, M-1 
DO 970 J=I+1,M 

ISUM2=ISUM2 + C MFC I )-l )*C MFC J )-l ) 

970 CONTINUE 
960 CC3NTINUE 

DO 965 1=1, M-2 

DO 975 J=I+1,M-1 
DO 985 K=J+1,M 

ISUM3=ISUM3 + C MFC I )-l )*C MFC J )-l )*C MFC K )-l ) 

985 CONTINUE 
975 CONTINUE 
965 CONTINUE 

ISUM=ISUM1+ISUM2+ISUM3 
IERDF=IT0T-ISUM-1 
MINR=MINIRC ITOT ,ISUM ) 

PRINT*, 'AT LEAST MINIMUM NUMBER OF REPLICATE PER TREATMENT =’,MINR 
PRINT*, ’DO YOU WANT TO ANALYSIS ONLY FACTOR-LEVEL C Y/N )? ’ 



59 



READ(5,7) ANS 

IF( (ANS.EQ. 'N' ).0R.( ANS.EQ. 'N* ) ) GO TO 990 

PRINT*, 'WHICH FACTOR ANALYSIS (IF YOU WANT 3 RO FACTOR, PRESS 3)?' 

READ(5,*) I 

NF1=MF(I)-1 

PICON=FLOAT( ITOT )/FLOAT( MF( I) ) 

GO TO 2000 

990 PRINT*, 'DO YOU WANT TO ANALYZE 2-WAY INTERACTION (Y/N)?' 

REA0(5,7) ANS 

IF( (ANS. E(3. 'N* ). OR. (ANS.EQ. 'N* )) GO TO 995 

PRINT*, 'WHICH INTERACTION ANALYSIS ?( IF YOU WANT 2 AND 3 WAY INTE 
&ACTION,THEN PRESS 2 3 ).' 

READ(5,*) I,J 

NF1 = (MF(I)-1)*(MF( J)-l) 

PICON=FLOAT( ITOT )/FLOAT( MF( I )*MF( J ) ) 

GO TO 2000 

995 PRINT*, 'DO YOU WANT TO ANALYZE 3-WAY INTERACTION (Y/N)?' 

READ(5,7) ANS 

IF( (ANS.EQ. 'N' ).OR.( ANS.EQ. 'N' ) ) GO TO 500 

PRINT*, 'WHICH INTERACTION ANALYSIS?( IF YOU WANT 2 AND 3 AND 4 INT 
&ERACTION ,THEN PRESS 2 3 4 ). ' 

READ(5,*) I, J, K 

NF1 = (MF( I )-l)*(MF( J)-l )*(MF(K)-1) 

PICON=FLOAT ( ITOT )/FLOAT( MF( I )*MF( J )*MF ( K ) ) 

GO TO 2000 

2000 PRINT*, 'DO YOU WANT TO PLOT R( ^ OF REPLICATE PER CELL) 

&VS POWER ( Y/N)?' 

READ(5,7) ANS 

IF( (ANS.EQ. 'Y' ). OR. (ANS.EQ. 'Y' )) GO TO 4100 

PRINT*, 'DO YOU WANT TO PLOT NONCENTRALITY VS POWER (Y/N)?' 

READ(5,7) ANS 

IF( (ANS.EQ. 'Y' ). OR. (ANS.EQ. 'Y' ) ) GO TO 4200 

PRINT*, 'DO YOU WANT TO PLOT ALPHA-LEVEL VS POWER (Y/N)?' 

READ(5,7) ANS 

IF(( ANS.EQ. 'Y' ).OR.( ANS.EQ. 'Y' ) ) GO TO 2300 
PRINT*, 'DO YOU WANT TO PLOT AGAIN (Y/N)?' 

READ (5,7) ANS 

IF( (ANS.EQ. 'Y' ). OR. (ANS.EQ. 'Y' )) GO TO 500 
GO TO 999 

C 

C 

C ANALYZE NUMBER OF SAMPLE SIZE PER TREATMENT VS POWER IN MULTI-WAY 
C ANALYSIS OF VARIANCE. 

C 

C 

4100 PRINT*, 'ALPHA LEVEL?' 

READ (5,*) ALPHA 

PRINT*, 'STANDARDIZED SQUARED DISTANCE VALUE ?' 

READ (5,*) STDS 

PRINT* ,' MAXIMUM R VALUE ( THE MAXIMUM VALUE ON X-AXIS )?' 

READ (5,*) MAXR 

PRINT*,' MINIMUM R VALUE ( MUST BE MORE THAN',MINR,' )?' 

READ (5,*) MINR 

PRINT*, 'INCREMENT R VALUE ?' 

READ (5,*) INCRR 
WRITE( 6,4095) 

4095 FORMAT( IX, 'DOFl DOF2 ALPHA-LEVEL F-INVERSE NONCENTRAL POWER') 
JJ = 1 

DO 4110 N=MINR, MAXR, INCRR 
NF2=N*ITOT-ISUM-l 
FINV=FINVER ( NFl ,NF2 , ALPHA ) 

NONCEN=N*PICON*STDS 

POW=FPO( ALPHA, NFl, NF2,FINV,NONCEN) 

WRITE ( 6 ,4015 ) NFl ,NF2 , ALPHA , FINV ,NONCEN ,POW 
4015 FORMAT( '0 ' ,I5,3X,I5,2X,2F10 .5,1X,F12. 2 ,1X,F10 .5 ) 

NN( JJ)=N 
POWER( JJ)=POW 
JJ=JJ+1 

4110 CONTINUE 
JJ=JJ-X 

PRINTit,'NN POWER' 

DO 4076 N=1,JJ 

WRITE( 6,4075 ) NN( N ) ,POWER( N ) 

4075 FORMAK 'O' , 'NN= ' , F8. 3 ,4X, ' POWER= ' , FIO .5 ) 

4076 CONTINUE 

CALL PLOTKNN, POWER, JJ,0) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN (Y/N)?' 

READ (5,7) ANS 

IF( (ANS.EQ. 'Y' ). OR. (ANS.EQ. 'Y' )) GO TO 500 
GO TO 999 

C 

C 

C ANALYZE NONCENTRALITY PARAMETER VS POWER IN MULTI-WAY 
C ANALYSIS OF VARIANCE. 

C 
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c 

4200 PRINT*, 'HOW MANY NUMNER OF REPLICATE (MUST BE MORE THANSMINR, 

)?• 

READ(5,*) IRR 
NF2=IRR*ITOT-ISUM-l 
PRINT* , • ALPHA-LEVE L ? ' 

READ (5,*) ALPHA 

PRINT*, 'MAXIMUM STANDARDIZED SQUARED DISTANCE VALUE RANGE ?' 

READ (5,*) MAXSD 

PRINT*, 'MINIMUM STANDARDIZED SQUARED DISTANCE VALUE RANGE ?' 

READ (5,*) MINSD 

PRINT*, 'INCREMENT STANDARDIZED SQUARED DISTANCE VALUE RANGE ?' 

READ (5,*) INCSD 
WRITE! 6,4495) 

4495 FORMATdX, ‘DOFl DOF2 ALPHA-LEVEL F-INVERSE NONCENTRAL POWER') 
JJ=1 

DO 4911 STDS=MINSD, MAXSD, INCSD 
FINV=FINVER ( NFl ,NF2 , ALPHA ) 

NONCEN=IRR*PICON*STDS 

POW=FPO( ALPHA, NFl, NF2,FINV,NONCEN) 

WRITE ( 6,4915 ) NFl ,NF2 , ALPHA , FINV,NONCEN ,POW 
4915 FORMAT( ' 0 ' ,I5,3X,I5,2X,2F10 .5,1X,F12. 2 ,1X,F10 .5 ) 

NN( JJ)=NONCEN 
POWER( JJ)=POW 
JJ=JJ+1 
4911 CONTINUE 
JJ=JJ-1 

DO 1476 I =1,JJ 

WRITE! 6,1475) NN( I ) ,POWER( I ) 

1475 FORMAT! 'O', 'NN=' ,F8.3,4X, 'POWER=' ,F10.5) 

1476 CONTINUE 

CALL PLOTT!NN, POWER, JJ,0) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN !Y/N)?' 

READ !5,7) ANS 

IF! ! ANS.EQ. 'Y' ).OR.! ANS.EQ. 'Y' )) GO TO 500 
GO TO 999 

C 

C 

C ANALYZE ALPHA-LEVEL PARAMETER VS POWER IN MULTI-WAY 
C ANALYSIS OF VARIANCE. 

C 

C 

2200 PRINT*, 'HOW MANY NUMNER OF REPLICATE (MUST BE MORE THAN' , MI NR, 

&• )?' 

READ! 5,*) IRR 
NF2=IRR*ITOT-ISUM-l 
READ !5,*) K 

PRINT*, 'STANDARDIZED SQUARED DISTANCE VALUE ?' 

READ (5,*) STDS 

PRINT*, 'MAXIMUM ALPHA RANGE ?' 

READ (5,*) MAXALP 

PRINT*, 'MINIMUM ALPHA RANGE ?' 

READ !5,*) MINALP 

PRINT*, 'INCREMENT ALPHA VALUE ?' 

READ (5,*) INCRAL 
WRITE (6,1295) 

1295 FORMAT! IX, 'DOFl DOF2 ALPHA-LEVEL F-iNVERSE NONCENTRAL POWER') 
JJ=1 
NF1=K-1 
NF2=K*!N-1) 

DO 1211 ALPHA=MINALP, MAXALP, INCRAL 
FINV=FINVER ! NFl ,NF2 , ALPHA ) 

NONCEN=IRR*PICON*STDS 

POW=FPO! ALPHA ,NF1 ,NF2 , FINV ,NONCEN ) 

WRITE! 6,1215) NFl, NF2, ALPHA, FINV, NONCEN,POW 
1215 FORMAT! '0 ' ,I5,3X,I5,2X,2F10.5,1X,F12.2,1X,F10.5) 

NN! JJ)=ALPHA 
POWER! JJ)=POW 
JJ=JJ+1 
1211 CONTINUE 

PRINT*> *NN( JJ ) ' ,MINK ,MAXK , INCRK 
DO 1276 I =1,JJ 

WRITE ! 6 , 12 75 ) NN! JJ ) , POWER ! J J ) 

1275 FORMAT! 'O', 'NN= ' ,F8. 3 ,4X, ' POWER= ' , FIO .5 ) 

1276 CONTINUE 

CALL PLOTT!NN, POWER, JJ,0) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN !Y/N)?' 

READ (5,7) ANS 

IF! (ANS.EQ. 'Y' ) .OR. ! ANS. EQ. 'Y' )) GO TO 500 
GO TO 999 



999 STOP 
END 
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* APPROXIMATION TO THE F« INVERSE DISTRIBUTION. * 

* CALCULATE F-INVERSE GIVEN DOFl ,D0F2 , ALPHA . * 

^ Q( F-ALPHA |D0F1,D0F2 )= ALPHA. * 

* * 



FUNCTION FINVER(N1,N2,ALP) 

REAL ALP,NUMER,DENUM,LAM,XP,YP,A,B,FIN1 
DATA CO, Cl, C2/2. 515517, 0.802853,0.010528/ 
DATA D1 ,D2,D5/1. -^52788, 0 . 189269,0 . 001308/ 
B=Nl/2. 

A=N2/2. 



CALCULATE NORMAL QUANTILE. 



T=SQRT(ALOG(l./(ALP*^^2. ))) 
NUMER=C0+C1*T+C2^^T*^^2 
DENUM=1 . +D1*T+D2*T**2+D3*T**3 
XP=T-(NUMER/DENUM) 



CALCULATE THE INVERSE F-DISTRIBUTION . 



LAM=(XP**2-3. >/6. 

W2= (l./(2.*B-l. ))-(!./( 2. *A-1. )) 

^l./(2.*A-l. )) + (!./( 2. *B-1. )) 

H = -1. ) 

Wl=( LAM+(5./6. )-( 2./(3.*H) ) ) 

W3=( XP*( H + LAM )**.5 )/H 
W=(W3-W2^^W1 )*2 

IF(W.LT.17^) THEN 
FIN1=EXP(W) 

ELSE 

FIN1=1000 

ENDIF 

FINVER=FIN1 

RETURN 

END 

GIVEN DOFl,DOF2,FINV,NON-CENTRAL PARAMETER, * 

THEN COMPUTE POWER. * 

( I.E. CALCULATE NON-CENTRAL F-DISTRIBUTION ( C.D.F.) 



FUNCTION FPO( ALP ,N1 ,N2 , FIN , NON ) 

REAL NON 
INTEGER N1,N2 

DATA A1 , A2, A3 , P/ . 4361836 , - . 1201676 , . 9372980 , . 33267/ 
CONST = 1.0 / SQRT (2.0 * 3.1415927) 



CALCULATE THE NORMAL APPROXIMATION . 



C1=(N1*FIN )/(Nl+NON) 

C2=l-( 2./( 9.*N2) ) 

C3=( 2,*(N1+2.*N0N) )/( 9.*(N1+N0N)**2) 

CCl = (Cl*^t(l./3. ))*C2-(1-C3) 

CC2 = ( (C3 + ( 2./( 9.^^N2) )*(Cl)**(2./3. ))**.5) 
FX=CC1/CC2 



CALCULATE THE NONCENTRAL F-DISTRIBUTION USING NORMAL CDF. 



CDF » 0.0 
IFCFX .EQ.O) THEN 
CDF = 0.5 

ELSE IF(FX.GT.O) THEN 
T=1./(1.+P*FX) 

XSQU=(FX^^FX)*.5 
IF(XSQU.LT.174) THEN 

2PDF=C0NST*EXP( -XSQU ) 

ELSE 

2PDF-0 

ENDIF 

C ZPDF=CONST*EXP( -XSQU ) 

CDF=1-2PDF*(A1*T+A2*T**2+A3^^T**5) 

ELSE 

FX=-FX 

T=1./(1.+P*FX) 
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XSQU=(FX*FX)*.5 
IF(XSQU.LT.174) THEN 

ZPDF=CONST*EXP( -XSQU ) 

ELSE 

ZPDF=0 

ENDIF 

CDF=ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

END IF 
FP0=1-CDF 
RETURN 
END 

•Jt •Jt 

* GIVEN ITOT AND ISUM , THEN CALCULATE MINIMUM REPLICATE. * 

* * 
**********************^(************************************************ 

FUNCTION MINIRl ITOT, ISUM) 

INTEGER IR, ITOT , ISUM, ICHECK,MINIR 
IR=2 

10 ICHECK=IR*IT0T-ISUM-1 
IF(ICHECK.LE.O) THEN 
IR=IR+1 
GO TO 10 
ELSE 

MINIR=IR 

ENDIF 

RETURN 

END 
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APPENDIX G 

THE POWER OF T-TEST PROGRAM LIST 



*********************************************************************** 



* * 

* PONER ALGORITHM ( T TEST ) * 

* * 

* DIRECTED BY : PROFESSOR DONALD. R. BARR * 

* * 

* WRITTEN BY : HUR, SEONG PIL JULY 1986 * 

* DEPARTMENT OF O.R. NAVAL POSTGRADUATE SCHOOL * 

* * 

* PROGRAM IS USED THE PRESENTATION OF POWER ALGORITHM OF T-TEST. * 

* (1) ONE-SAMPLE CASE TEST. * 

* . ONE-SIDED TEST ( HA : MU>MU0 ) * 

* . ONE-SIDED TEST ( HA : MU<MUO ) * 

* . TWO-SIDED TEST ( HA : MU MUD ) * 

* (2) TWO-SAMPLE CASE TEST. * 

* . ONE-SIDED TEST ( HA : MU>MU0 ) * 

* . ONE-SIDED TEST ( HA : MU<MU0 ) * 

* . TWO-SIDED TEST ( HA : MU MUD ) * 

* NOTE : Subroutine PLOTT is NONIMSL subroutine library. * 

* * 



REAL NN( 200 ) ,POWER( 200 ) ,DF ,NONCDF ,TCDF1 ,TCDF2 ,TINV1 
REAL TINV,NONCEN,TX, ALPHA, CDF, CONST ,MAXMU1 ,MINMU1 ,INCMU1 
REAL MAXDEL,MINDEL,INCDEL,P0W,P0W1,P0W2,MUZER0,STDV,MU1 
INTEGER JJ,K,N 
CHARACTER*! ANS 

C 



c 

c 

c 

c 


INPUT USER OPTIONS. 




100 


PRINT*, 'DO YOU WANT TO TEST ONE SAMPLE 
READ(5,10) ANS 


T-TEST(YZN)?' 


10 


FORMAT( Al) 






IF( (ANS.EQ. *N' ). OR. (ANS.EQ. 'N' ) ) GO TO 
PRINT*,' SAMPLE SIZE ?' 

READ(5,*) N 

PRINT*,' ALPHA-LEVEL ?' 

READ(5,*) ALPHA 
PRINT*,' MU-ZERO VALUE ?' 

READ(5,*) MUZERO 

PRINT*,' STANDARD DEVIATION VALUE ?' 


160 


c 


READ(5,*) STDV 





c 

C ONE -SAMPLE ONE-SIDED TEST, HA : MU > MUO. 

C 

C 

PRINT*, 'DO YOU WANT TO USE ONE-SIDED TEST (Y/N)?' 

RE AD (5,10) ANS 

IF( (ANS.EQ. 'N* ).OR.(ANS.EQ. 'N* )) GO TO 50 

PRINT*, 'DO YOU WANT TO TEST H0:MU=MU0 VS HA:MU>MU0 (Y/N)?' 

READ(5,10) ANS 

IF( (ANS.EQ. 'N* ). OR. (ANS. E(3. *N' ) ) GO TO ^0 

PRINT*,' MAXIMUM MUl VALUE ( THE MAXIMUM VALUE OF X-AXIS )?' 
READ (5,*) MAXMUl 

PRINUt,* MINIMUM MUl VALUE ( THE MINIMUM VALUE OF X-AXIS )? 
a (CONDITION ;MU1 MUST BE MORE THAN MU-2ER0 )' 

READ (5,*) MINMUl 

PRINTif, • LARGE INCREMENT N VALUE ?' 

READ (5,«) INCMUl 
N=N-1 

C 

c 

C GENERATE T-INVERSE VALUE ( CRITICAL REGION VALUE ). 

C 

C 

TINV=TINVER( N , ALPHA ) 

JJ=1 
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DO 50 MU1=MINMU1,MAXMU1,INCMU1 



C 

C 

C COMPUTE NON-CENTRALITY VALUE. 

C 

C 

DF=N 

NONCEN=( (MU1-MUZER0)*SQRT(DF+1 ) )/STDV 

C 

C 

C GENERATE CDF OF NONCENTRAL T DISTRIBUTION. 

C 

C 

TCDF=NONCDF(N,NONCEN,TINV) 

POW=l-TCDF 
NN( JJ)=MU1 
POWER( JJ)=POW 
JJ=JJ+1 

50 CONTINUE 

JJ=JJ-1 

PRINT*, 'NN POWER' 

DO 60 1=1, JJ 

WRITE (6,70 ) NN(I),POWER(I) 

70 FORMAT! 'O' , 'NN= ' , F8 . 3 ,^X, ' POWER= ' ,F10 . 5 ) 

60 CONTINUE 

CALL PLOTTINN, POWER, JJ,0) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN (Y/N)?' 

READ (5,10) ANS 

IF((ANS.EQ. 'Y* ).OR.(ANS.EQ. 'Y' )) GO TO 100 
GO TO 80 

C 

c 

C ONE-SAMPLE ONE-SIDED TEST, HA : MU < MUO. 

C 

C 

^0 PRINT*, 'DO YOU WANT TO TEST HO:MU=MUO VS HA:MU<MUO (Y/N)?' 
READ(5,10) ANS 

IF( ( ANS.EQ. 'N* ).0R.( ANS.EQ. 'N' ) ) GO TO 80 

PRINT*,' MAXIMUM MUl VALUE ( THE MAXIMUM VALUE OF X-AXIS )? ' 
READ (5,*) MAXMUl 

PRINT*,' MINIMUM MUl VALUE ( THE MINIMUM VALUE OF X-AXIS )? 

& (CONDITION :MU1 MUST BE LESS THAN MU-ZERO ) ' 

READ (5,*) MINMUl 

PRINT*,' INCREMENT MUl VALUE ?' 

READ (5,*)INCMU1 
N=N-1 

TINV=TINVER( N, ALPHA ) 

TINV=-TINV 

JJ=1 

DO 90 MU1=MINMU1, MAXMUl, INCMUl 
DF=N 

NONCEN=( ( MUl-MUZERO )*SQRT( DF + 1 ) )/STDV 
TCDF=NONCDF(N,NONCEN,TINV) 

POW=TCDF 
NN( JJ)=MU1 
POWER! JJ)=POW 
JJ=JJ+1 
90 CONTINUE 
JJ=JJ-1 

PRINT*, 'NN POWER* 

DO 110 1=1, JJ 

WRITE (6,120) NN( I ) , POWER! I ) 

120 FORMAT! '0*, 'NN= * ,F8.3,^X, 'POWER= * ,F10.S) 

110 CONTINUE 

CALL PLOTKNN, POWER, JJ,0) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN (Y/N)?* 

READ (5,10) ANS 

IF! (ANS.EQ. 'Y* ). OR. (ANS.EQ. *Y' )) GO TO 100 
GO TO 80 

C 

c 

C ONE-SAMPLE TWO-SIDED TEST, HA : MU =/ MUO. 

C 

C 

30 PRINT*, 'DO YOU WANT TO USE TWO-SIDED TEST! Y/N)?* 

READ«5,10) ANS 

IF((ANS.EQ. 'N* ). OR. (ANS.EQ. 'N' ) ) GO TO 80 

PRINT*, 'MAXIMUM MUl VALUE ( THE MAXIMUM VALUE ON X-AXIS )?* 
READ (5,*) MAXMUl 

PRINT*,* MINIMUM MUl VALUE ( THE MINIMUM VALUE ON X-AXIS )? 

& (CONDITION :MU1 MUST BE LESS THAN MU-ZERO )* 

READ (5,*) MINMUl 

PRINT*,' INCREMENT MUl VALUE ?* 

READ (5,*)INCMU1 
N=N-1 
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ALPHA=ALPHA/2. 

TINV=TINVER( N , ALPHA ) 

JJ=1 

DO 130 MU1=MINMU1,MAXMU1,INCMU1 
DF=N 

NONCEN=nnUl-MUZERO)^SQRT(DF + l. ) )/STDV 
TCDFl=NONCDF(N,NONCEN,TINV) 

P0W1=1.-TCDF1 

TINV1=-TINV 

TCDF2=NONCDF( N ,NONCEN ,TINT71 ) 

POH2=TCDF2 
POW=POWl+POW2 
NN( JJ)=MU1 
POWER( JJ )=POW 
JJ=JJ+1 
130 CONTINUE 
JJ=JJ-1 

PRINT*, *NN POWER' 

DO 1^0 1=1, JJ 

WRITE (6,150) NN(I ),POWER(I ) 

150 FORMAT! 'O' , 'NN=* ,F8.3,4X, 'POWER= ' ,F10 .5 ) 

1^0 CONTINUE 

CALL PLOTTINN, POWER, JJ, 0 ) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN (Y/N)?' 

READ (5,10) ANS 

IF( ( ANS.EQ. 'Y* ).OR. (ANS.EQ. 'Y' ) ) GO TO 100 
GO TO 80 

C 

C 

C TWO-SAMPLE ONE-SIDED TEST, HA : MUX > MUY. 

C 

C 

160 PRINT*, 'DO YOU WANT TO TEST TWO SAMPLE T-TEST( Y/N )? ' 

READ! 5,10) ANS 

IF! ( ANS.EQ. 'N' ).OR.(ANS.EQ. *N' )) GO TO 80 
PRINT*,' SAMPLE SIZE FROM POPULATION 1 ?' 

READ!5,*) N 

PRINT*,' SAMPLE SIZE FROM POPULATION 2 ?' 

READ!5,*) M 

PRINT*,' ALPHA-LEVEL ?' 

READ! 5,*) ALPHA 

PRINT*, 'DO YOU WANT TO USE ONE-SIDED TEST !Y/N)?' 

READ (5,10) ANS 

IF!!ANS.EQ. 'N' ). OR. (ANS.EQ. 'N' )) GOTO 170 

PRINT*, 'DO YOU WANT TO TEST HO: MUX=MUY VS HA: MUX>MUY (Y/N)?' 
READ! 5,10) ANS 

IF! (ANS.EQ. 'N' ). OR. (ANS.EQ. 'N' )) GO TO 180 
PRINT*,' MAXIMUM NONCENTRALITY VALUE ? ' 

READ (5,*) MAXDEL 

PRINT*,' MINIMUM NONCENTRALITY VALUE?' 

C & (CONDITION :MU1 MUST BE MORE THAN MU-ZERO )? ' 

READ (5,*) MINDEL' 

PRINT*,' INCREMENT NONCENTRALITY VALUE?' 

READ (5,*)INCDEL 
Nl=M+N-2 

TINV=TINVER( N1 , ALPHA ) 

JJ = 1 

DO 190 NONCEN=MINDEL, MAXDEL, INCDEL 
TCDF=NONCDF( N1 ,NONCEN ,TINV ) 

POW=l-TCDF 
NN( JJ)=NONCEN 
POWER! JJ)=POW 
JJ=JJ+1 
190 CONTINUE 
JJ=JJ-1 

PRINT*, 'NN POWER* 

DO 200 1=1, JJ 

WRITE(6,210) NN(I),POWER(I ) 

210 FORMAT! 'O' > 'NN= ' , F8 . 3 ,4X, * POWER= ' ,F10 . 5 ) 

200 CONTINUE 

CALL PLOTT( NN, POWER, JJ,0) 

PRINT^f>'DO YOU WANT TO PLOT AGAIN (Y/N)?' 

READ (5,10) ANS 

IF! (ANS.EQ. 'Y' ). OR. (ANS.EQ. 'Y' )) GO TO 100 
GO TO 80 

C 

C 

C TWO-SAMPLE ONE-SIDED TEST, HA : MUX < MUY. 

C 

C 

180 PRINT*, 'DO YOU WANT TO TEST H0;MUX=mY VS HA:MUX<MUY (Y/N)?' 
READ! 5,10) ANS 

IF((ANS.EQ. 'N' ). OR. (ANS.EQ. 'N' )) GO TO 80 
PRINT*,' MAXIMUM NONCENTRALITY VALUE ? ' 

READ (5,*) MAXDEL 
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PRINT*,' MINIMUM NONCENTRALITY VALUE?' 

REAO (5,*) MINDEL 

PRINT*,' INCREMENT NONCENTRALITY VALUE?' 

REAO (5,*)INCOEL 
Nl=M+N-2 

TINV=TINVER( N1 , ALPHA ) 

TINV2=-TINV 

JJ=1 

DO 220 NONCEN=MINDEL,MAXDEL,INCDEL 
TCDF=NONCDF(Nl,NONCEN,TINV2) 

POH=TCOF 
NN( JJ)=NONCEN 
POWER( JJ)=POW 
JJ=JJ+1 

220 CONTINUE 
JJ=JJ-1 

PRINT*, 'NN POHER' 

DO 230 1=1, JJ 

WRITE( 6,240) NN(I ),POWER(I ) 

240 FORMAT('0', 'NN= ‘ ,F8. 3 ,4X , ' POWER= ' ,F10 . 5 ) 

230 CONTINUE 

CALL PLOTKNN, POHER, JJ,0) 

PRINT*, 'DO YOU WANT TO PLOT AGAIN (Y/N)?' 

READ (5,10) ANS 

IF((ANS.EQ.'Y' ).OR.(ANS.EQ.'Y' )) GO TO 100 
GO TO 80 

C 

C 

C TWO-SAMPLE TWO-SIOEO TEST, HA ; MUX =/ MUY . 

C 

C 

170 PRINT*, 'DO YOU WANT TO USE TWO-SIOED TEST (Y/N)?' 
READ(5,10) ANS 

IF((ANS.EQ. 'N' ).OR.(ANS.EQ. 'N' )) GO TO 80 
PRINT*,’ MAXIMUM NONCENTRALITY VALUE ? ' 

READ (5,*) MAXDEL 

PRINT*,' MINIMUM NONCENTRALITY VALUE?' 

REAO (5,*) MINDEL 

PRINT* , • INCREMENT NONCENTRALITY VALUE? ' 

READ (5,*)INC0EL 

Nl=M+N-2 

ALPHAl=ALPHA/2. 

TINV=TINVER( N1 , ALPHAl ) 

JJ=1 

00 250 NONCEN=MINOEL,MAXOEL,INCOEL 
TCDF1=N0NCDF(N1,N0NCEN,TINV) 

P0W1=1.-TCDF1 

TINV1=-TINV 

TCDF2=NONCOF(N1,NONCEN,TINV1) 

POW2=TCDF2 
POW=POWl+POW2 
NN( JJ)=NONCEN 
POWER( JJ)=POW 
JJ=JJ+1 
250 CONTINUE 
JJ=JJ-1 

PRINT*, 'NN POWER' 

00 260 1 = 1, JJ 

WRITE( 6,270) NN( I ) ,P0WER( I ) 

270 FORMAT('0', 'NN= ' ,F8.3,4X, 'POWER=' ,F10.5) 

260 CONTINUE 

CALL PLOTT(NN,POWER,JJ,0) 

PRINT*, '00 YOU WANT TO PLOT AGAIN (Y/N)?' 

REAO (5,10) ANS 

IF((ANS.EQ. 'Y' ).OR.(ANS.EQ. 'Y' )) 60 TO 100 
GO TO 80 
80 STOP 
ENO 



* * 

* APPROXIMATION TO THE T- INVERSE OISTRIBUTION. * 

* CALCULATE T-INVERSE GIVEN OOF ALPHA. * 

* Q( T-ALPHA I OOF )= ALPHA. * 

* 



FUNCTION TINVER(N, ALPHA) 

REAL ALPHA, NUMER,OENUM, LAM, XP,G1,G2, 63, X,TINVER 
INTEGER N 

DATA CO, Cl, C2/2. 515517, 0.802853, 0.010328/ 

OATA 01,02,03/1.432788,0.189269,0.001308/ 

C 

C 

C CALCULATE NORMAL QUANTILE. 

C 

C 
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T=SQRT( ALOG( 1 ./( ALPHA**2. ) ) ) 

NUMEfx=C0+Cl*T+C2*T**2 

DENUM=1 . +D1*T+D2*T**2 . +D3*T**3 . 

XP=T-(NUMER/DENUM) 

C 

C 

C CALCULATE THE INVERSE T-DISTRIBUTION. 

C 

C 

X=XP 

G1=.25*(X**3+X) 

G2 = ( 1 . /96 . )*( 5 . *X**5+16 . *X**3+3 . *X ) 

G3 = ( 1 . /ZS ^ . )* ( 3 . *X** 7+ 1 9 . *X**3 + 17. *X**3 - 15 . *X ) 

G^= ( 1 . /92160 . )* ( 79 . *X**9+776 . *X**7+1482 . *X**5-1592 . *X**3-945 . *X ) 

T P =X + ( G 1/N ) + ( G 2/ ( N*N ) ) ♦ ( G3/ ( N^* 3 ) ) + ( G^/ ( ) ) 

TINVER=TP 

RETURN 

END 

*********************** INFORMATION ********************************** 



* * 

* THE OBJECTIVE OF THIS PROGRAM IS TO CALCULATE NON-CENTRAL T * 

* DISTRIBUTION (CDF) * 

* THE CALCULATION METHOD IS NORMAL APPROXIMATION. * 

* * 



********************** VARIABLE DECLARATION ************************** 

FUNCTION NONCDF(NF,TLAM,TSTAR) 

INTEGER NF 

REAL TSTAR,TLAM,TX, NONCDF,CDF 

DATA A1 , A2 , A3 , P/ . 4361836 , - . 1201676 , . 9372980 , . 33267/ 

CONST = 1.0 / SORT (2.0 * 3.1415927) 

C 

C 

C CALCULATE THE NORMAL APPROXIMATION . 

C 

C 

Cl=( !.-(!./( 4. *NF) ) )*TSTAR-TLAM 
C2=SQRT( 1 . +( TSTAR**2/( 2 .*NF ) ) ) 

TX=C1/C2 
CDF = 0.0 
IF(TX .EQ.O) THEN 
CDF = 0.5 

ELSE IF(TX.GT.O) THEN 
T=1./(1.+P*TX) 

XSQU=(TX*TX)*.5 
IF(XSQU.LT.174) THEN 

ZPDF=CONST*EXP( -XSQU ) 

ELSE 

ZPDF=0 

ENDIF 

CDF=1-ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

ELSE 

TX=-TX 

T=1./(1.+P*TX) 

XSQU=(TX*TX)*.5 
IF( XSQU. LT. 174) THEN 

ZPDF=CONST*EXP( -XSQU) 

ELSE 

ZPDF=0 

ENDIF 

CDF=ZPDF*( A1*T+A2*T**2+A3*T**3 ) 

END IF 

NONCDF=CDF 

RETURN 

END 
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